Student Member, IEEE Jie Chen

Member, IEEE Paul Honeine

Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model

Spectral unmixing is an important issue to analyze remotely sensed hyperspectral data. Although the linear mixture model has obvious practical advantages, there are many situations in which it may not be appropriate and could be advantageously replaced by a nonlinear one. In this paper, we formulate a new kernel-based paradigm that relies on the assumption that the mixing mechanism can be described by a linear mixture of endmember spectra, with additive nonlinear fluctuations defined in a reproducing kernel Hilbert space. This family of models has clear interpretation, and allows to take complex interactions of endmembers into account. Extensive experiment results, with both synthetic and real images, illustrate the generality and effectiveness of this scheme compared with state-of-the-art methods.

I. INTRODUCTION

Hyperspectral imaging is a continuously growing area of remote sensing, which has received considerable attention in the last decade. Hyperspectral data provide a wide spectral range, coupled with a high spectral resolution. These characteristics are suitable for detection and classification of surfaces and chemical elements in the observed images. Applications include land use analysis, pollution monitoring, wide-area reconnaissance, and field surveillance, to cite a few. Due to multiple factors, including the possible low spatial resolution of some hyperspectral-imaging devices, the diversity of materials in the observed scene, the reflections of photons onto several objects, etc., mixed-pixel problems can occur and be critical for proper interpretation of images. Indeed, assigning mixed pixels to a single pure component, or endmember, inevitably leads to a loss of information.

Spectral unmixing is an important issue to analyze remotely sensed hyperspectral data. This involves the decomposition of each mixed pixel into its pure endmember spectra, and the estimation of the abundance value for each endmember [START_REF] Keshava | Spectral unmixing[END_REF]. Several approaches have been developed for endmember extraction [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. On the one hand, methods with pure pixel assumption have been proposed to extract the endmembers from pixels in the scene, such as the pixel purity index algorithm [START_REF] Boardman | Atomatic spectral unmixing of AVIRIS data using convex geometry concepts[END_REF], the vertex component analysis (VCA) [START_REF] Nascimento | Vertex Component Analysis: A fast algorithm to unmix hyperspectral data[END_REF], and the N-FINDR algorithm [START_REF] Winter | N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF], among others [START_REF] Plaza | Recent developments in endmember extraction and spectral unmixing[END_REF], [START_REF] Honeine | Geometric unmixing of large hyperspectral images: A barycentric coordinate approach[END_REF]. On the other hand, some methods have been proposed to overcome the absence of pure pixels, by generating virtual endmembers, such as the minimum volume simplex analysis (MVSA) [START_REF] Li | Minimum volume simplex analysis: A fast algorithm to unmix hyperspectral data[END_REF], the minimum volume enclosing simplex algorithm (MVES) [START_REF] Chan | A convex analysisbased minimum-volume enclosing simplex algorithm for hyperspectral unmixing[END_REF], and the minimum volume constrained nonnegative matrix factorization (MVC-NMF) [START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[END_REF]. Endmember identification and abundance estimation can be conducted either in a sequential or collaborative manner. Under the assumption that the endmembers have been identified, hyperspectral image unmixing then reduces to estimating the fractional abundances. The term unmixing in the paper represents the abundance estimation step, which is referred to as the supervised unmixing in some literature.

The linear mixture model is widely used to identify and quantify pure components in remotely sensed images due to its simple physical interpretation and trackable estimation process. To be physically interpretable, the driving abundances are often required to satisfy two constraints: all abundances must be nonnegative, and their sum must be equal to one. In addition to the extremely low-complexity method that has been recently proposed [START_REF] Honeine | Geometric unmixing of large hyperspectral images: A barycentric coordinate approach[END_REF], which is based on geometric considerations, at least two classes of approaches can be distinguished to determine abundances. On the one hand, there are estimation methods that lead to an optimization problem which must be solved subject to non-negativity and sum-to-one constraints [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF]. On the other hand, following the principles of Bayesian inference, there are simulation techniques that define prior distributions for abundances, and estimate unknown parameters based on the resulting joint posterior distribution [START_REF] Moussaoui | On the decomposition of Mars hyperspectral data by ICA and bayesian positive source separation[END_REF], [START_REF] Dobigeon | Joint bayesian endmember extraction and linear unmixing for hyperspectral imagery[END_REF], [START_REF] Eches | Bayesian estimation of linear mixtures using the normal compositional model. application to hyperspectral imagery[END_REF], [START_REF] Themelis | A novel hierarchical bayesian approach for sparse semisupervised hyperspectral unmixing[END_REF]. Some recent works also take sparsity constraints into account in the unmixing process [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF], [START_REF] Themelis | A novel hierarchical bayesian approach for sparse semisupervised hyperspectral unmixing[END_REF], [START_REF] Guo | L1 unmixing and its application to hyperspectral image enhancement[END_REF], [START_REF] Bioucas-Dias | Hyperspectral unmixing: Geometrical, statistical, and sparse regression-based approaches[END_REF], [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF].

Although the linear mixture model has obvious practical advantages, there are many situations in which it may not be appropriate (e.g., involving multiple light scattering effects) and could be advantageously replaced by a nonlinear one. For instance, multiple scattering effects can be observed on complex vegetated surfaces [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] where it is assumed that incident solar radiation is scattered by the scene through multiple bounces involving several endmembers. Some nonlinear mixture models, such as the generalized bilinear model studied in [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], account for presence of multi-photon interactions by introducing additional interaction terms in the linear model. Another typical situation is the case where the components of interest are in an intimate association, and the photons interact with all the materials simultaneously as they are multiply scattered. A bidirectional reflectance model based on the fundamental principles of radiative transfer theory was proposed in [START_REF] Hapke | Bidirectional reflectance spectroscopy, 1, Theory[END_REF] to describe these interactions. It is usually referred to as the intimate mixture model. Obviously, the mixture mechanism in a real scene may be much more complex than the above models and often relies on scene parameters that are difficult to obtain.

Nonlinear unmixing has generated considerable interest among researchers, and different methods have been pro-posed to account for nonlinear effects. Using training-based approaches is a way to bypass difficulties with unknown mixing mechanism and parameters. In [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF], a radial basis function neural network was used to unmix intimate mixtures. In [START_REF] Plaza | Nonlinear neural network mixture models for fractional abundance estimation in AVIRIS hyperspectral images[END_REF], the authors designed a multi-layer perceptron neural network combined with a Hopfield neural network to deal with nonlinear mixtures. In [START_REF] Plaza | Recent developments in endmember extraction and spectral unmixing[END_REF], [START_REF] Plaza | On the use of small training sets for neural network-based characterization of mixed pixels in remotely sensed hyperspectral images[END_REF], the authors discussed methods for automatic selection and labeling of training samples. These methods require the networks to be trained using pixels with known abundances, and the quality of the training data may affect the performance notably. Moreover, for a new set of spectra in a scene, or different embedded parameters, a new neural network should be trained again before unmixing can be performed. Approaches that do not require training samples were also studied in the literature. In [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], a nonlinear unmixing algorithm for the general bilinear mixture model was proposed. Based on Bayesian inference, this method however has a high computational complexity and is dedicated to the bilinear model. In [START_REF] Raksuntorn | Nonlinear spectral mixture analysis for hyperspectral imagery in an unknown environment[END_REF], [START_REF] Nascimento | Nonlinear mixture model for hyperspectral unmixing[END_REF], the authors extended the collection of endmembers by adding artificial cross-terms of pure signatures to model light scattering effects on different materials. However, it is not easy to identify which cross-terms should be selected and added to the endmember dictionary. If all the possible cross-terms were considered, the set of endmembers would expand dramatically. In [START_REF] Nascimento | Unmixing hyperspectral intimate mixtures[END_REF], the authors addressed the nonlinear unmixing problem with an intimate mixture model. The proposed method first converts observed reflectances into albedo using a look-up table, then a linear algorithm estimates the endmember albedos and the mass fractions for each sample. This method is based on the hypothesis that all the parameters of the intimate mixture model are known. Nonlinear algorithms operating in reproducing kernel Hilbert spaces (RKHS) have been a topic of considerable interest in the machine learning community, and have proved their efficiency in solving nonlinear problems. Kernel-based methods have already been considered for detection and classification in hyperspectral images [START_REF] Camps-Valls | Kernel-based methods for hyperspectral image classification[END_REF], [START_REF] Heesung | Kernel orthogonal subspace projection for hyperspectral signal classification[END_REF]. Kernel-based nonlinear unmixing approaches have also been investigated [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF], [START_REF] Broadwater | A comparison of kernel functions for intimate mixture models[END_REF], [START_REF] Wu | A kernel spatial complexity-based nonlinear unmixing method of hyperspectral imagery[END_REF]. These algorithms were mainly obtained by replacing each inner product between endmember spectra, in the cost functions to be optimized, by a kernel function. This can be viewed as a nonlinear distortion map applied to the spectral signature of each material, independently of their interactions. This principle may be extremely efficient in solving detection and classification problems as a proper distortion can increase the detectability or separability of some patterns. It is however of little physical interest in solving the unmixing problem because the nonlinear nature of mixing is not only governed by individual spectral distortions, but also by nonlinear interactions of the materials.

In this paper, we formulate the problem of estimating abundances of a nonlinear mixture of hyperspectral data. This new kernel-based paradigm allows to take nonlinear interactions of the endmembers into account. It leads to a more meaningful interpretation of the unmixing mechanism than existing kernel based methods. The abundances are determined by solving an appropriate kernel-based regression problem under constraints. This paper is organized as follows. Section II introduces the basic concepts of our modeling approach. Section III presents a new kernel-based hyperspectral mixture model, called K-Hype, and the associated identification algorithm to extract the abundances within this nonlinear context. The balance between linear and nonlinear contributions is unfortunately fixed in K-Hype. In order to overcome this drawback, a natural generalization called SuperK-Hype or SK-Hype, is then largely described. It is based on the concept of Multiple Kernel Learning, and allows to automatically adapt the balance between linear spectral interactions and nonlinear ones. Finally, major differences with some existing works on kernelbased processing of hyperspectral images are also pointed out. In Section IV, experiments are conducted using both synthetic and real images. Performance comparisons with other popular methods are also provided. Section V concludes this paper and gives a short outlook onto future work.

II. A KERNEL-BASED NONLINEAR UNMIXING PARADIGM

Let r = [r 1 , r 2 , . . . , r L ] ⊤ be an observed column pixel, supposed to be a mixture of R endmember spectra, with L the number of spectral bands. Assume that M = [m 1 , m 2 , . . . , m R ] is the (L × R) target endmember matrix, where each column m i is an endmember spectral signature. For the sake of convenience, we shall denote by m ⊤ λ ℓ the ℓ-th (1 × R) row of M , that is, the vector of the R endmember signatures at the ℓ-th wavelength band. Let α = [α 1 , α 2 , . . . , α R ] ⊤ be the abundance column vector associated with the pixel r.

We first consider the linear mixing model where any observed pixel is a linear combination of the endmembers, weighted by the fractional abundances, that is,

r = M α + n ( 1
)
where n is a noise vector. Under the assumption that the endmember matrix M is known, the vector α of fractional abundances is usually determined by minimizing a cost function of the form

J(α) = J reg (α) + 1 2µ r -M α 2 = J reg (α) + 1 2µ L ℓ=1 (r ℓ -α ⊤ m λ ℓ ) 2 (2)
under the non-negativity and sum-to-one constraints1

α i ≥ 0, ∀i ∈ 1, . . . , R R i=1 α i = 1, (3)
where J reg (•) is a regularization function, and µ is a small positive parameter that controls the trade-off between regularization and fitting. The above analysis assumes that the relationship between m λ ℓ and r ℓ is dominated by a linear function. There are however many situations, involving multiple scattering effects, in which model (1) may be inappropriate and could be advantageously replaced by a nonlinear one.

Consider the general mixing mechanism

r = Ψ(M ) + n (4)
with Ψ an unknown nonlinear function that defines the interactions between the endmembers in matrix M . This requires us to consider a more general problem of the form

ψ * = arg min ψ∈H 1 2 ψ 2 H + 1 2µ L ℓ=1 (r ℓ -ψ(m λ ℓ )) 2 (5)
with H a given functional space, and µ a positive parameter that controls the trade-off between regularity of the function ψ and fitting. Clearly, this basic strategy may fail if the functionals ψ of H cannot be adequately and finitely parameterized.

Kernel-based methods rely on mapping data from the original input space into a feature space by means of a nonlinear function, and then solving a linear problem in that new space. They lead to efficient and accurate resolution of the inverse problem [START_REF] Winter | N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF], as it has been showed in the literature. See, e.g., [START_REF] Schölkopf | Advances in kernel methods[END_REF], [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. Our paper exploits the central idea of this research area, known as the kernel trick, to investigate new nonlinear unmixing algorithms. We shall now review the main definitions and properties related to reproducing kernel Hilbert spaces [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] and Mercer kernels [START_REF] Mercer | Functions of positive and negative type and their connection with the theory of integral equations[END_REF].

Let H denote a Hilbert space of real-valued functions ψ on a compact M, and let • , • H be the inner product in the space H. Suppose that the evaluation functional δ m λ defined by δ m λ [ψ] = ψ(m λ ) is linear with respect to ψ and bounded, for all m λ in M. By virtue of the Riesz representation theorem, there exists a unique positive definite function m λ → κ(m λ , m λj ) in H, denoted by κ(•, m λj ) and called representer of evaluation at m λj , which satisfies [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] ψ(m λj ) = ψ, κ(•, m λj ) H , ∀ψ ∈ H [START_REF] Plaza | Recent developments in endmember extraction and spectral unmixing[END_REF] for every fixed m λj ∈ M. A proof of this may be found in [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. Replacing ψ by κ(•, m λi ) in ( 6) yields

κ(m λi , m λj ) = κ(•, m λi ), κ(•, m λj ) H (7)
for all m λi , m λj ∈ M. Equation ( 7) is the origin of the generic term reproducing kernel to refer to κ. Denoting by Φ the map that assigns the kernel function κ(•, m λj ) to each input data m λj , equation [START_REF] Honeine | Geometric unmixing of large hyperspectral images: A barycentric coordinate approach[END_REF] immediately implies that κ(m λi , m λj ) = Φ(m λi ), Φ(m λj ) H . The kernel thus evaluates the inner product of any pair of elements of M mapped to the space H without any explicit knowledge of Φ and H. Within the machine learning area, this key idea is known as the kernel trick.

The kernel trick has been widely used to transform linear algorithms expressed only in terms of inner products into nonlinear ones. Considering again [START_REF] Winter | N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF], the optimum function ψ * can be obtained by solving the following least squares support vector machines (LS-SVM) problem [START_REF] Suykens | Least Squares Support Vector Machines[END_REF]

ψ * = arg min ψ∈H 1 2 ψ 2 H + 1 2µ L ℓ=1 e 2 ℓ subject to e ℓ = r ℓ -ψ(m λ ℓ ), ℓ ∈ {1, . . . , L} (8)
We introduce the Lagrange multipliers β ℓ , and consider the Lagrange function

L = 1 2 ψ 2 H + 1 2µ L ℓ=1 e 2 ℓ - L ℓ=1 β ℓ (e ℓ -r ℓ + ψ(m λ ℓ )). (9)
The conditions for optimality with respect to the primal variables are given by

ψ * = L ℓ=1 β * ℓ κ(•, m λ ℓ ) e * ℓ = µ β * ℓ ( 10
)
We then derive the dual optimization problem

β * = arg max β - 1 2 β ⊤ (K + µI) β + β ⊤ r, ( 11
)
where K is the so-called Gram matrix whose (ℓ, p)-th entry is defined by κ(m λ ℓ , m λp ). Classic examples of kernels are the radially Gaussian kernel

κ(m λ ℓ , m λp ) = exp -m λ ℓ -m λp 2 /2σ 2 , and the Laplacian kernel κ(m λ ℓ , m λp ) = exp -m λ ℓ -m λp /σ , with σ ≥ 0 the kernel bandwidth. Another example of interest is the q-th degree polynomial kernel κ(m λ ℓ , m λp ) = (1 + m ⊤ λ ℓ m λp ) q , with q ∈ IN * .
The kernel function κ maps m λ ℓ into a very high, even infinite, dimensional space H without any explicit knowledge of the associated nonlinear function. The vector β * and κ then describe the relation between the endmembers and the observation. The goal of the analysis is however to estimate the abundance vector, and there is no direct relation between α * and β * in the general case. In what follows, we shall focus attention on the design of specific kernels that enable us to determine abundance fractions within this context.

III. KERNEL DESIGN AND UNMIXING ALGORITHMS

The aim of this section is to propose kernel design methods and the corresponding algorithms to estimate abundances. The two approaches described hereafter are flexible enough to capture wide classes of nonlinear relationships, and to reliably interpret a variety of experimental measurements. Both have clear interpretation.

A. A preliminary approach for kernel-based hyperspectral unmixing: the K-Hype algorithm

In order to extract the mixing ratios of the endmembers, we define the function ψ in (5) by a linear trend parameterized by the abundance vector α, combined with a nonlinear fluctuation term, namely,

ψ(m λ ℓ ) = α ⊤ m λ ℓ + ψ nlin (m λ ℓ ) subject to α 0 and 1 ⊤ α = 1 (12)
where ψ nlin can be any real-valued functions on a compact M, of a reproducing kernel Hilbert space denoted by H nlin . Let κ nlin be its reproducing kernel. It can be shown [START_REF] Haussler | Convolution kernels on discrete structures[END_REF] that, as the direct sum H lin ⊕ H nlin of the RKHS of kernels

κ lin (m λ ℓ , m λp ) = m ⊤ λ ℓ m λp and κ nlin (m λ ℓ , m λp ) defined on
M, the space H of functions of the form ( 12) is also a RKHS with kernel function

κ(m λ ℓ , m λp ) = (κ lin ⊕ κ nlin )(m λ ℓ , m λp ) = m ⊤ λ ℓ m λp + κ nlin (m λ ℓ , m λp ). (13)
The corresponding Gram matrix K is given by

K = M M ⊤ + K nlin ( 14
)
where K nlin is the Gram matrix associated with the nonlinear map ψ nlin , with (ℓ, p)-th entry κ nlin (m λ ℓ , m λp ).

We propose to conduct hyperspectral data unmixing by solving the following convex optimization problem

ψ * = arg min ψ 1 2 ψ lin 2 Hlin + ψ nlin 2 Hnlin + 1 2µ L ℓ=1 e 2 ℓ where ψ = ψ lin + ψ nlin with ψ lin (m λ ℓ ) = α ⊤ m λ ℓ subject to e ℓ = r ℓ -ψ(m λ ℓ ) α 0 and 1 ⊤ α = 1 (15)
By the strong duality property, we can derive a dual problem that has the same solution as the above primal problem. Let us introduce the Lagrange multipliers β ℓ , γ r and λ. The Lagrange function associated with the problem (15) can be written as

G = 1 2 α 2 + ψ nlin 2 Hnlin + 1 2µ L ℓ=1 e 2 ℓ - L ℓ=1 β ℓ (e ℓ -r ℓ + ψ(m λ ℓ )) - R r=1 γ r α r + λ(1 ⊤ α -1) (16)
with γ r ≥ 0. We have used that ψ lin 2 Hlin = α 2 because the functional space H lin , parametrized by α, contains all the function of the variable m λ ℓ of the form ψ lin (m λ ℓ ) = α ⊤ m λ ℓ . It is characterized by the norm

ψ lin 2 Hlin = κ lin (α, • ), κ lin (α, • ) H = κ lin (α, α) = α 2 .
(17) The conditions for optimality of G with respect to the primal variables are given by

   α * = L ℓ=1 β * ℓ m λ ℓ + γ * -λ * 1 ψ * nlin = L ℓ=1 β * ℓ κ nlin (•, m λ ℓ ) e * ℓ = µ β * ℓ ( 18
)
By substituting ( 18) into ( 16), we get the dual problem [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] (next page). Provided that the coefficient vector β * has been determined, the measured pixel can be reconstructed using

r * = [ψ * (m λ1 ), . . . , ψ * (m λL )] ⊤ = M (M ⊤ β * + γ * -λ * 1) + K nlin β * (20)
as indicated by [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF]. Comparing the above expression with [START_REF] Moussaoui | On the decomposition of Mars hyperspectral data by ICA and bayesian positive source separation[END_REF], we observe that the first and the second term of the r.h.s. of equation ( 20) correspond to the linear trend and the nonlinear fluctuations, respectively. Finally, the abundance vector α * can be estimated as follows

α * = M ⊤ β * + γ * -λ * 1 (21)
Problem ( 19) is a quadratic program (QP). Numerous candidate methods exist to solve it, such as interior point, active set and projected gradient, as presented in [START_REF] Bertsekas | Nonlinear programming[END_REF], [START_REF] Luenberger | Linear and nonlinear programming[END_REF]. These well known numerical procedures lie beyond the scope of this paper.

B. Some remarks on kernel selection

Selecting an appropriate kernel is of primary importance as it captures the nonlinearity of the mixture model. Though an infinite variety of possible kernels exists, it is always desirable to select a kernel that is closely related to the application context. The following example justifies the combination [START_REF] Moussaoui | On the decomposition of Mars hyperspectral data by ICA and bayesian positive source separation[END_REF], which associates a linear model with a nonlinear fluctuation term. It also allows us to define a possible family of appropriate kernels for data unmixing.

Consider the generalized bilinear mixing model presented in [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], at first, limited to three endmember spectra for the sake of clarity

r =M α + γ 12 α 1 α 2 (m 1 ⊗ m 2 ) + γ 13 α 1 α 3 (m 1 ⊗ m 3 ) + γ 23 α 2 α 3 (m 2 ⊗ m 3 ) + n with α 1 , α 2 , α 3 ≥ 0 α 1 + α 2 + α 3 = 1 (22)
where γ 12 , γ 13 and γ 23 are attenuation parameters, and ⊗ the Hadamard product. It can be observed that the nonlinear term with respect to α, in the r.h.s. of [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF], is closely related to the homogeneous polynomial kernel of degree 2, that is,

κ nlin (m λ ℓ , m λp ) = (m ⊤ λ ℓ m λp ) 2 .
Indeed, with a slight abuse of notation, the latter can be written in an inner product form as follows

κ nlin (m λ ℓ , m λp ) = Φ nlin (m λ ℓ ) ⊤ Φ nlin (m λp ) (23)
with

Φ nlin (m λ ℓ ) =(m 2 λ ℓ,1 , m 2 λ ℓ,2 , m 2 λ ℓ,3 , √ 2m λ ℓ,1 m λv,2 , √ 2m λ ℓ,1 m λ ℓ,3 , √ 2m λ ℓ,2 m λ ℓ,3 ) ⊤ (24)
where m λ ℓ,i is the i-th entry of m λ ℓ . This means that, in addition to the linear mixture term M α, the auto and interaction terms considered by the kernel-based model are of the form m i ⊗ m j for all i, j = 1, . . . , R.

By virtue of the reproducing kernel machinery, endmember spectra do not need to be explicitly mapped into the feature space. This allows to consider complex interaction mechanisms by changing the kernel κ nlin , without having to modify the optimization algorithm described in the previous subsection. As an illustration, consider the polynomial kernel

κ nlin (m λ ℓ , m λp ) = (1 + m ⊤ λ ℓ m λp ) q .
Making use of the binomial theorem yields

κ nlin (m λ ℓ , m λp ) = q k=0 q k (m ⊤ λ ℓ m λp ) k . ( 25
)
We observe that each component

(m ⊤ λ ℓ m λp ) k = (m λ ℓ,1 m λp,1 + . . . + m λ ℓ,R m λp,R ) k of the above expression max β,γ,λ G ′ (β, γ, λ) = - 1 2   β γ λ   ⊤   K + µI M -M 1 M ⊤ I -1 -1 ⊤ M ⊤ -1 ⊤ R     β γ λ   +   r 0 -1   ⊤   β γ λ   subject to γ 0 (19)
with

K = M M ⊤ + K nlin .
can be expanded into a weighted sum of k-th degree monomials of the form

(m λ ℓ,1 m λp,1 ) k1 (m λ ℓ,2 m λp,2 ) k2 . . . (m λ ℓ,R m λp,R ) kR (26) with R r=1 k r = k.
This means that, in addition to the linear mixture term M α, the auto and interaction terms considered by the kernel-based model are of the form m k1 1 ⊗ m k2 2 ⊗ . . . ⊗ m kR R for every set of exponents in the Hadamard sense satisfying 0 ≤ R r=1 k r ≤ q. Note that it would be computationally expensive to explicitly form these interaction terms. Their number is indeed very large: there are p k monomials (26) of degree k, and then 1-R q+1 1-R components in the entire q-th order representation. Compared with the methods introduced in [START_REF] Raksuntorn | Nonlinear spectral mixture analysis for hyperspectral imagery in an unknown environment[END_REF], [START_REF] Nascimento | Nonlinear mixture model for hyperspectral unmixing[END_REF], which insert products of pure material signatures as new endmembers, we do not need to extend the endmember matrix by adding such terms. The kernel trick makes the computation much more tractable. In the experimentations reported hereafter, the following 2-nd degree polynomial kernel was used

κ nlin (m λ ℓ , m λp ) = 1 + 1 R 2 (m λ ℓ -1/2) ⊤ (m λp -1/2) 2 ( 27
)
where the constants 1/R 2 and 1/2 serve the purpose of normalization.

C. Nonlinear unmixing by multiple kernel learning: the SK-Hype algorithm

The proposed model relies on the assumption that the mixing mechanism can be described by a linear mixture of endmember spectra, with additive nonlinear fluctuations ψ nlin defined in a RKHS. This justifies the use of a Gram matrix of the form K = M M ⊤ + K nlin in the algorithm presented previously. Model [START_REF] Moussaoui | On the decomposition of Mars hyperspectral data by ICA and bayesian positive source separation[END_REF] however has some limitations in that the balance between the linear component α ⊤ m λ ℓ and the nonlinear component ψ nlin (m λ ℓ ) cannot be tuned. This should however be made possible as recommended by physicallyinspired models such as model [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF]. In addition, kernels κ nlin with embedded linear component such as the inhomogeneous polynomial kernel (25) introduces a bias into the estimation of α, unless correctly estimated and removed. Another difficulty is that the model ( 12) cannot captures the dynamic of the mixture, which requires that r or the m λ ℓ 's be locally normalized. This unlikely situation occurs, e.g., if a library of reflectance signatures is used for the unmixing process. In order to address problems such as the above, we shall now consider Gram matrices of the form

K u = uM M ⊤ + (1 -u)K nlin ( 28
)
with u in [0, 1] in order to ensure positiveness of K u .

Learning both the parameter u and the mixing coefficients β in a single optimization problem is known as the multiple kernel learning problem. See [START_REF] Rakotomamonjy | SimpleMKL[END_REF] and references therein. The rest of this section is devoted to the formulation and resolution of a convex optimization problem.

In order to tune the balance between ψ lin (m λ ℓ ) = α ⊤ m λ ℓ and ψ nlin (m λ ℓ ), it might seem tempting to substitute matrix K with matrix K u in the dual problem [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF]. Unfortunately, a primal problem must be first formulated in order to identify, in the spirit of equation ( 18), explicit expressions for ψ lin and ψ nlin . We propose to conduct hyperspectral data unmixing by solving the following primal problem

ψ * , u * = arg min ψ,u 1 2 1 u ψ lin 2 Hlin + 1 1 -u ψ nlin 2 Hnlin + 1 2µ L ℓ=1 e 2 ℓ
subject to e ℓ = r ℓ -ψ(m λ ℓ ) and 0 ≤ u ≤ 1 where ψ = ψ lin + ψ nlin [START_REF] Heesung | Kernel orthogonal subspace projection for hyperspectral signal classification[END_REF] where u allows to adjust the balance between ψ lin and ψ nlin via their norms. The spaces H ′ lin and H ′ nlin are RKHS of the general form

H ′ = ψ ∈ H : ψ H u < ∞ (30)
with the convention x 0 = 0 if x = 0, and ∞ otherwise. This implies that, if u = 0, then ψ belongs to space H ′ if and only if ψ = 0. By continuity consideration via this convention, it can be shown that problem [START_REF] Heesung | Kernel orthogonal subspace projection for hyperspectral signal classification[END_REF] is a convex optimization problem by virtue of the convexity of the so-called perspective function f (u, ψ) = ψ 2 H ′ /u over IR + × H ′ . This has been shown in [START_REF] Boyd | Convex Optimization[END_REF]Chapter 3] in the finite-dimensional case, and extended in [START_REF] Rakotomamonjy | SimpleMKL[END_REF] to the infinite-dimensional case. This allows to formulate the two-stage optimization procedure, with respect to ψ and u successively, in order to solve problem [START_REF] Heesung | Kernel orthogonal subspace projection for hyperspectral signal classification[END_REF].

min u J(u) subject to 0 ≤ u ≤ 1 (31)
where J(u) is defined by (32) (next page). The connection between [START_REF] Heesung | Kernel orthogonal subspace projection for hyperspectral signal classification[END_REF] and this problem is as follows. We always have [42, p. 133]

min u,ψ F (u, ψ) = min u J(u) (33)
where J(u) = min ψ F (u, ψ), subject to all the constraints over u and ψ defined in ( 31)-( 32). In addition, as proven in

J(u) =            min ψ F (u, ψ) = 1 2 1 u ψ lin 2 Hlin + 1 1 -u ψ nlin 2 Hnlin + 1 2µ L ℓ=1 e 2 ℓ subject to e ℓ = r ℓ -ψ(m λ ℓ ) with ψ = ψ lin + ψ nlin and ψ lin (m λ ℓ ) = h ⊤ m λ ℓ with h 0 (32)
textbooks [42, p. 87], because F is convex in (u, ψ) subject to convex constraints over ψ, it turns out that J(u) is convex in u and, as a consequence, that the constrained optimization problem ( 31) is convex.

Compared to the preliminary algorithm described in Section III-A, it is important to note that the sum-to-one constraint 1 ⊤ h = 1 has been given up. We have adopted an alternative strategy, which consists of relaxing this constraint, and normalizing h afterward in order to be correctly interpreted as a vector α of fractional abundances. The main reasons are as follows. The weight vector h is optimally scaled by the factor u in order to achieve the best balance between ψ lin and ψ nlin . As a consequence, there is no significant advantage of satisfying the sum-to-one constraint 1 ⊤ h = 1 along the optimization process, which is consequently simplified. On the contrary, we clearly observed that relaxing the ℓ 1 -norm 1 ⊤ h of the weight vector h acts as an additional degree of freedom to adjust the balance between ψ lin and ψ nlin , which would otherwise be limited to the factors u and 1 -u. This has been confirmed by experiments, which have revealed a significant improvement in performance.

By the strong duality property, we shall now derive a dual problem that has the same solution J(u) = F (u, ψ * ) as the primal problem [START_REF] Wu | A kernel spatial complexity-based nonlinear unmixing method of hyperspectral imagery[END_REF]. Let us introduce the Lagrange multipliers β ℓ and γ r . The Lagrange function associated with the problem (32) can be written as

G = 1 2 1 u h 2 + 1 1 -u ψ nlin 2 Hnlin + 1 2µ L ℓ=1 e 2 ℓ - L ℓ=1 β ℓ (e ℓ -r ℓ + ψ(m λ ℓ )) - R r=1 γ r h r (34)
with γ r ≥ 0, where we have used that ψ lin 2 Hlin = h 2 . The conditions for optimality of G with respect to the primal variables are given by

     h * = u L ℓ=1 β * ℓ m λ ℓ + γ * ψ * nlin = (1 -u) L ℓ=1 β * ℓ κ nlin (•, m λ ℓ ) e * ℓ = µ β * ℓ ( 35
)
By substituting ( 35) into (34), we get the dual problem (36) (next page). Pixel reconstruction can be performed using

r * = [ψ * (m λ1 ), . . . , ψ * (m λL )] ⊤ with ψ * (m λ ℓ ) = m ⊤ λ ℓ h * + ψ * nlin (m λ ℓ
) defined in equation [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. Finally, the estimated abundance vector is given by

α * = M ⊤ β * + γ * 1 ⊤ (M ⊤ β * + γ * ) (37)
Let us briefly address the differentiability issue of the problem ( 31)- [START_REF] Mercer | Functions of positive and negative type and their connection with the theory of integral equations[END_REF]. The existence and computation of the derivatives of supremum functions such as J(u) have been largely discussed in the literature. As pointed out in [START_REF] Rakotomamonjy | SimpleMKL[END_REF], [START_REF] Bonnans | Optimization problems with perturbations: A guided tour[END_REF], the differentiability of J at any point u 0 is ensured by the unicity of the corresponding minimizer (β * 0 , γ * 0 ), and by the differentiability of the cost function F (u, ψ) in [START_REF] Wu | A kernel spatial complexity-based nonlinear unmixing method of hyperspectral imagery[END_REF]. The derivative of J at u 0 can be calculated as if the minimizer (β * 0 , γ * 0 ) was independent of u 0 , namely, dJ(u) du u=u0 = I summarizes the proposed algorithm. Note that ( 31) is a very small-size problem. Indeed, it involves a one-dimension optimization variable and can thus be solved with an ad-hoc procedure. Using a gradient projection method, e.g., based on Armijo rule along the feasible direction, makes practical sense in this case [START_REF] Bertsekas | Nonlinear programming[END_REF]Chapter 2]. See also [START_REF] Chen | Nonnegative least-mean-square algorithm[END_REF]. Moreover, both problems can benefit of warm-starting between successive solutions to speed-up the optimization procedure. The algorithm can be stopped based on conditions for optimality in convex optimization framework. In particular, the KKT conditions and the duality gap should be equal to zero, within a numerical error tolerance specified by the user. The variation of the cost J(u) between two successive iterations should also be considered as a potential stopping criterion.

∂G ′ (u,β * 0 ,γ * 0 ) ∂u u=u0 . This yields dJ(u) du u=u0 = - 1 2 M ⊤ β * 0 + γ * 0 2 -β * 0 ⊤ K nlin β * 0 (38) Table
Before testing our algorithms, and comparing their performance with state-of-the-art approaches, we shall now explain how they differ from existing kernel-based techniques for hyperspectral data processing.

D. Comparison with existing kernel-based methods in hyperspectral imagery

Some kernel-based methods have already been proposed to process hyperspectral images, with application to classification, supervised or unsupervised unmixing, etc. By taking advantage of capturing nonlinear data dependences, some of them have been shown to achieve better performance than their linear counterpart. Let us now briefly discuss the main difference between our kernel-based model and those presently existing. The central idea underlying most of state-of-theart methods is to nonlinearly transform hyperspectral pixelvectors prior to applying a linear algorithm, simply by replacing inner products with kernels in the cost function. This basic principle is fully justified in detection/classification problems because a proper nonlinear distortion of spectral signatures can increase the detectability/separability of materials. Within

J(u) =          max β,γ G ′ (u, β, γ) = - 1 2 β γ ⊤ K u + µI uM uM ⊤ uI β γ + r 0 ⊤ β γ subject to γ 0 with K u = uM M ⊤ + (1 -u)K nlin ( 36
)

Initialization

Choose the kernel κnlin and the regularization constant µ ≥ 0.

Calculate the kernel matrix Knlin.

Repeat

Calculate J(u) in (36) by using a generic QP solver Solve the optimization problem minu J(u) subject to 0 ≤ u ≤ 1 by performing one iteration of the gradient projection algorithm, using [START_REF] Haussler | Convolution kernels on discrete structures[END_REF] until stopping criterion is satisfied Estimate the abundances by [START_REF] Suykens | Least Squares Support Vector Machines[END_REF] the context of hyperspectral unmixing, this leads to consider mixtures of the form

Φ(r) = [Φ(m 1 ), Φ(m 2 ), . . . , Φ(m R )] α + n. ( 39
)
This model is inherent in the KFCLS algorithm [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF], [START_REF] Broadwater | A comparison of kernel functions for intimate mixture models[END_REF], which optimizes the following mean-square error criterion where all the inner products have been replaced by kernels

J(α) = α ⊤ K m α -2 α ⊤ κ rm + κ(r, r), (40)
where K m is the Gram matrix with (i, j)-th entry κ(m i , m j ), and κ rm is a vector with i-th entry κ(r, m i ). Unfortunately, even though model [START_REF] Bertsekas | Nonlinear programming[END_REF] allows distortions of spectral signatures, it does not explicitly include nonlinear interactions of the endmember spectra. The analysis in Section III-B has shown strong connections between our kernel-based model and wellcharacterized models, e.g., the generalized bilinear mixture model. The experimental comparison on simulated and real data reported in the next section confirms this view.

IV. EXPERIMENTAL RESULTS

We shall now conduct some simulations to validate the proposed unmixing algorithms, and to compare them with state-of-the-art methods, using both synthetic and real images.

A. Experiments on synthetic images

Let us first report some experimental results on synthetic images, generated by linear and nonlinear mixing of several endmember signatures. The materials we have considered are alunite, calcite, epidote, kaolinite, buddingtonite, almandine, jarosite and lepidolite. They were selected from the ENVI software library. These spectra consist of 420 contiguous bands, covering wavelengths ranging from 0.3951 to 2.56 micrometers.

In the first scene, only three materials were selected to generate images: epidote, kaolinite, buddingtonite. In the second scene, five materials were used: alunite, calcite, epidote, kaolinite, buddingtonite. In the third scene, the eight materials were used. For each scene, three 50-by-50 hyperspectral images were generated with different mixture models, each providing N = 2500 pixels for evaluating and comparing the performance of several algorithms. These three models were the linear model ( 1), the bilinear mixture model defined as

r = M α + R-1 i=1 R j=i+1 α i α j (m i ⊗ m j ) + n, (41)
and a post-nonlinear mixing model (PNMM) [START_REF] Jutten | Advances in nonlinear blind source separation[END_REF] defined by

r = (M α) ξ + n (42)
where ( • ) ξ denotes the exponential value ξ applied to each entry of the input vector. Parameter ξ was set to 0.7. The abundance vectors α n , with n = 1, . . . , 2500, were uniformly generated in the simplex defined by non-negative and sum-toone constraints. Finally, all these images were corrupted with an additive white Gaussian noise n with two levels of SNR, 30 dB and of 15 dB. The following algorithms were considered

• The so-called Fully Constrained Least Square method (FCLS), [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF]: This technique was derived based on linear mixture model. It provides the optimal solution in the least-mean-square sense, subject to non-negativity and sum-to-one constraints. A relaxation parameter ν has to be tuned to specify a compromise between the residual error and the sum-to-one constraint. • The extended endmember-matrix method (ExtM), [START_REF] Raksuntorn | Nonlinear spectral mixture analysis for hyperspectral imagery in an unknown environment[END_REF]: This method consists of extending the endmember matrix M artificially with cross-spectra of pure materials in order to model light scatter effects. In the experiments, all the second-order cross terms m i ⊗ m j were inserted so that it would correspond to the generalized bilinear model. This approach also has a relaxation parameter ν for the sum-to-one constraint.

• The so-called Kernel Fully Constrained Least Square method (KFCLS), [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF]: This is a kernel method, directly derived from FCLS, in which all the inner products are replaced by kernel functions. As for all the other kernelbased algorithms considered in this paper, the Gaussian kernel was used for simulations. This algorithm has two parameters, the bandwidth σ of the Gaussian kernel, and a relaxation parameter ν for the sum-to-one constraint. • The Bayesian algorithm derived for generalized bilinear model (BilBay), [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF]: This method is based on appropriate prior distributions for the unknown parameters, which must satisfy the non-negativity and sum-to-one constraints, and then derives joint posterior distribution of these parameters. A Metropolis-within-Gibbs algorithm is used to estimate the unknown model parameters. The MMSE estimates of the abundances were computed by averaging the 2500 generated samples obtained after 500 burn-in iterations.

• The first algorithm proposed in this paper (K-Hype):

This is the preliminary algorithm described in Section III-A. The Gaussian kernel (G) with bandwidth σ, and the polynomial kernel (P) defined by [START_REF] Nascimento | Unmixing hyperspectral intimate mixtures[END_REF] were considered. The Matlab optimization function Quadprog was used to solve the QP problem. • The second algorithm proposed in this paper (SK-Hype): This is the main algorithm described in Section III-C and Table I. As for K-Hype, the Gaussian kernel and the polynomial kernel were considered. In order to simplify the experiments, the weight vector u was constrained to be of the form u = [u, 1 -u] ⊤ with 0 ≤ u ≤ 1. Obviously, this allows to reduce the number of variables but does not change the structure of the algorithm. The variable u was initially set to 1 2 . A gradient projection method, based on the Armijo rule to compute the optimal step size along the feasible direction, was used to determined u. The algorithm was stopped when the relative variation of u between two successive iterations became less than ζ max = 10 -3 , or the maximum number of iterations It max = 10 was reached. The Matlab optimization function Quadprog was used to solve the QP problem. The root mean square error defined by

RMSE = 1 N R N n=1 α n -α * n 2 (43)
was used to compare these six algorithms. In order to tune their parameters, preliminary runs were performed on 100 independent test pixels for each experiment. The bandwidth σ of the Gaussian kernel in the algorithms ExtM, K-Hype and SK-Hype was varied within {1, . . . , 3} with increment of 1 2 . The parameter µ of K-Hype and SK-Hype algorithms was varied within {1, 10 -1 , 10 -2 , 5 • 10 -3 }. The parameter ν in algorithms FCLS, ExtM, KFCLS was chosen within {1, 10 -1 , 10 -2 , 10 -3 , 10 -4 }. All the parameters used in the experiments are reported in the Appendix.

Results for Scene 1 to Scene 3 unmixing, with three, five and eight endmember materials, are reported in Table II, Table III and Table IV respectively. Because the FCLS method was initially derived for the linear mixing model, it achieves a very low RMSE for linearly-mixed images, and produces a relatively large RMSE with nonlinearly-mixed images. With second-order cross terms that extend the endmember matrix M , the ExtM algorithm notably reduces the RMSE when dealing with bilinearly-mixed images when compared with FCLS. However, it marginally improves the performance in PNMM image ummixing. BilBay algorithm was derived for the bilinear mixing model, and thus achieves very good performance with bilinearly-mixed images. Nevertheless, the performance of BilBay clearly degrades when dealing with a nonlinear mixing model for which it was not originally designed. KFCLS with Gaussian kernel performs worse than FCLS, even with nonlinearly-mixed images as it does not clearly investigate nonlinear interactions between materials.

For the less noisy scenes (30 dB), our algorithms K-Hype and SK-Hype exhibit significantly reduced RMSE when dealing with nonlinearly-mixed images. In the case of the bilinear model, K-Hype and SK-Hype achieve very good performance compared to the other algorithms. Indeed, they are the best performers except in a few cases. In the case of the PNMM model, they outperform all the other algorithms, and it can be observed that SK-Hype outperforms K-Hype in several scenarios. For the noisiest scenes (15 dB), although the increase in the noise level significantly degrades the performance of all the algorithms, K-Hype and SK-Hype still maintain an advantage. Last but not least, the margin of performance over the other approaches becomes larger as the number of endmembers increases.

To give a more meaningful comparison of the performance of these algorithms, one-tailed Welch's t-tests with significance level 0.05 were used to test the hypothesis

H 0 : RMSE proposed < RMSE literature
where RMSE proposed denotes the RMSE of the K-Hype and SK-Hype algorithms, with Gaussian and polynomial kernels, and RMSE literature is the RMSE of the algorithms of the literature selected in this paper. Due to limited space, only the results for Scene 2 and the SNR level 30 dB are reported here, in Table V to VII. The letter A means that the hypothesis H 0 is accepted. Without ambiguity, these results confirm the advantage of our algorithms.

FCLS ExtM KFCLS BilBay K-Hype (G)

A K-Hype (P) A SK-Hype (G) A SK-Hype (P) A

The computational time of these algorithms mainly depends on the constrained optimization problem to be solved. FCLS and KFLCS minimize a quadratic cost function of dimensionR, under inequality constraints of the same dimension. ExtM solves a similar problem but with an increased

FCLS ExtM KFCLS BilBay K-Hype (G)

A A A A K-Hype (P) A A A SK-Hype (G) A A A A SK-Hype (P) A A A A TABLE VII WELSH'S t-TESTS FOR SCENE 2 WITH SNR = 30 DB (PNMM) FCLS ExtM KFCLS BilBay K-Hype (G) A A A A K-Hype (P) A A A A SK-Hype (G) A A A A SK-Hype (P) A A A A
dimension due to the cross-spectra that are artificially inserted.

In the case where only the second-order cross spectra are added, the dimension of the optimization problem is R + R 2 with R = 3, 5 and 8 in this study. BilBay has to generate numerous samples to estimate the model parameters, and suffers from the large computational cost of this sampling strategy. K-Hype solves a quadratic programming problem of dimension L + R + 1. It is interesting to note that the computational cost is independent of the complexity of the unmixing model. A sparsification strategy as described in [START_REF] Richard | Online prediction of time series data with kernels[END_REF] should be advantageously used to greatly reduce the computational complexity with negligible effect on the quality of the results. SK-Hype has similar advantages as K-Hype except that the alternating optimization scheme requires more time. The average computational times per pixel of all these algorithms are listed in Table VIII. 2

B. Experiment with AVIRIS image

This section illustrates the performance of the proposed algorithms, and several other algorithms, when applied to real hyperspectral data. The scene that was used for our experiment is the well-known image captured on the Cuprite mining district (NV, USA) by AVIRIS. A sub-image of 250 × 191 pixels was chosen to evaluate the algorithms. This area of interest has L = 188 spectral bands. The number of endmembers was first estimated via the virtual dimensionality, and R was accordingly set to 12 [START_REF] Nascimento | Vertex Component Analysis: A fast algorithm to unmix hyperspectral data[END_REF]. VCA algorithm was then used to 2 Using Matlab R2008a on a iMac with 3.06GHz Intel Core 2 Duo and 4 Go Memory extract the endmembers. Both our algorithms were compared with all the state-of-the-art algorithms considered previously.

L = 420 L = 210 R = 3 R = 5 R = 8 R = 3 R = 5 R = 8 FCLS 7.
After preliminary experiments, the regularization parameters of FCLS and ExtM algorithms were set to ν = 0.01. K-Hype algorithm and SK-Hype algorithm were run with the polynomial kernel [START_REF] Nascimento | Unmixing hyperspectral intimate mixtures[END_REF], and the Gaussian kernel. The bandwidth of the Gaussian kernel was set to σ = 2. The regularization parameter was fixed to 2 • 10 -3 . To evaluate the performance, the averaged spectral angle between original r and reconstructed r * pixel vectors was used

Θ = 1 N N n=1 θ(r n , r * n )
where N is the number of processed pixels and θ(r, r * ) = cos -1 r,r * r r *

. It is important to note that the quality of reconstruction, estimated by the averaged spectral angle or mean-square error for instance, is not necessarily in proportion to the the quality of unmixing, especially for real images where the nonlinear mixing mechanism can be complex. In particular, more complicated model may better fit the data. Parameter Θ is only reported here as complementary information. The averaged spectral angle of each approach is reported in Table IX. Note that KFCLS was not considered in these tests as there is no possible direct reconstruction of pixels. Clearly, our algorithms have much lower reconstruction errors than the other approaches. Six typical estimated abundance maps out of twelve available are shown in Figure 1. It can be observed that the estimated locations of the different materials are quite similar for the four methods, except the US Highway 95 in the last column which is much more accurately depicted by our methods. Finally, the distributions of reconstruction errors r nr * n /L associated to these methods are shown in Figure2. V. CONCLUSION Spectral unmixing is an important issue to analyze remotely sensed hyperspectral data. This involves the decomposition of each mixed pixel into its pure endmember spectra, and the estimation of the abundance value for each endmember. To be physically interpretable, the abundances are often required to be nonnegative, and their sum must be equal to one. Although the linear mixture model has many practical advantages, there are many situations in which it may not be appropriate. In this paper, we formulated a new kernel-based paradigm that relies on the assumption that the mixing mechanism can be described by a linear mixture of endmember spectra, with additive nonlinear fluctuations defined in a RKHS. This family of models has a clear physical interpretation, and allows to take complex interactions of endmembers into account. Two kernel-based algorithms for estimating the abundances were proposed. The second one, based on the concept of multiple kernel learning, is the natural generalization of the first one as it allows to automatically adapt the balance between linear spectral interactions and nonlinear ones. Future work will include studying the feasibility and constraints of designing physically meaningful kernels, possibly based on manifold learning as in [START_REF] Heylen | Non-linear spectral unmixing by geodesic simplex volume maximization[END_REF], [START_REF] Nguyen | Hyperspectral image unmixing using manifold learning methods: derivations and comparative tests[END_REF]. We will also focus our attention on adaptive kernel-based algorithms, in the spirit of [START_REF] Richard | Online prediction of time series data with kernels[END_REF], to unmix neighboring pixel-vectors sequentially and thus speed up processing.

Fig. 1 .

 1 Fig. 1. Abundances maps of selected materials. From top to bottom: FCLS, BilBay, K-Hype (G), SK-Hype (G). From left to right: calcedony, alunite, kaolinite, buddingtonite, sphene, US highway 95.

Fig. 2 .

 2 Fig. 2. Maps of reconstruction error. From left to right: FCLS, BilBay, K-Hype (G), SK-Hype (G).

TABLE I NONLINEAR

 I UNMIXING BY MULTIPLE KERNEL LEARNING: THE SK-HYPE ALGORITHM

TABLE V WELSH

 V 'S t-TESTS FOR SCENE 2 WITH SNR = 30 DB (LINEAR MODEL)

TABLE VI WELSH

 VI 'S t-TESTS FOR SCENE 2 WITH SNR = 30 DB (BILINEAR MODEL)

TABLE VIII AVERAGED

 VIII COMPUTATIONAL TIME PER PIXEL (IN SECONDS)

  6 × 10 -4 1.2 × 10 -3 2.4 × 10 -3 6.4 × 10 -4 1.1 × 10 -3 2.1 × 10 -3 ExtM 1.6 × 10 -3 5.5 × 10 -3 0.020 1.4 × 10 -3 4.0 × 10 -3 0.014 KFCLS 6.0 × 10 -4 1.3 × 10 -3 2.5 × 10 -3 5.6 × 10 -4 8.4 × 10 -4 1 × 10 -3

	BilBay	6.6	15	40	6.4	14	33
	K-Hype	0.17	0.29	0.47	0.038	0.065	0.11
	SK-Hype	4.3	7.5	13.4	0.92	1.7	3.1

J. Chen and P. Honeine are with the Institut Charles Delaunay, CNRS, Université de Technologie de Troyes, France (e-mail: chenjieg@sina.com, paul.honeine@utt.fr) C. Richard is with the Université de Nice Sophia-Antipolis, CNRS, Observatoire de la Côte d'Azur, France (e-mail: cedric.richard@unice.fr)

For ease of notation, these two constraints will be denoted by α 0 and 1 ⊤ α = 1, where 1 is a vector of ones.