Fei Zhu
email: fei.zhu@utt.fr

Abderrahim Halimi
email: a.halimi@hw.ac.uk

Paul Honeine
email: paul.honeine@univ-rouen.fr

Badong Chen
email: chenbd@mail.xjtu.edu.cn

Nanning Zheng
email: nnzheng@mail.xjtu.edu.cn

ADMM for Maximum Correntropy Criterion

Keywords: Correntropy, maximum correntropy estimation, alternating direction method of multipliers, hyperspectral image, unmixing problem

The correntropy provides a robust criterion for outlier-insensitive machine learning, and its maximisation has been increasingly investigated in signal and image processing. In this paper, we investigate the problem of unmixing hyperspectral images, namely decomposing each pixel/spectrum of a given image as a linear combination of other pixels/spectra called endmembers. The coefficients of the combination need to be estimated subject to the nonnegativity and the sum-to-one constraints. In practice, some spectral bands suffer from low signal-to-noise ratio due to acquisition noise and atmospheric effects, thus requiring robust techniques for the unmixing problem. In this work, we cast the unmixing problem as the maximization of a correntropy criterion, and provide a relevant solution using the alternating direction method of multipliers (ADMM) method. Finally, the relevance of the proposed approach is validated on synthetic and real hyperspectral images, demonstrating that the correntropybased unmixing is robust to outlier bands.

I. INTRODUCTION

Spectral unmixing is an essential issue in many disciplines, including machine learning and signal/image processing, with a wide range of applications, such as classification, segmentation, material identification and target detection. Typically, a hyperspectral image corresponds to a scene taken at many continuous and narrow bands across a certain wavelength range; namely, each pixel is a spectrum. Assuming that each spectrum is a mixture of several pure materials, the unmixing problem consists in identifying these pure materials (the socalled endmembers) and estimating their proportions (the so-called abundances) at each pixel [START_REF] Keshava | Spectral unmixing[END_REF]. While the endmember extraction is relatively easy from geometry (see [START_REF] Honeine | Geometric unmixing of large hyperspectral images: a barycentric coordinate approach[END_REF] and references therein), the abundance estimation is still an open problem. Indeed, the abundances can be estimated using leastsquares methods, geometric approaches [START_REF]A simple scheme for unmixing hyperspectral data based on the geometry of the n-dimensional simplex[END_REF], [START_REF] Honeine | Geometric unmixing of large hyperspectral images: a barycentric coordinate approach[END_REF], or by including nonlinearities with kernel-based machines [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF], [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF], [START_REF] Nguyen | Supervised nonlinear unmixing of hyperspectral images using a preimage method[END_REF], [START_REF] Chen | Nonlinear estimation of material abundances of hyperspectral images with 1 -norm spatial regularization[END_REF], [START_REF]Estimating abundance fractions of materials in hyperspectral images by fitting a post-nonlinear mixing model[END_REF], [START_REF] Chen | Nonlinear unmixing of hyperspectral data with partially linear least-squares support vector regression[END_REF]. In this paper, we consider the abundance estimation in unmixing hyperspectral images.

The linear mixture model (LMM) is the most investigated over the past decades [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF], [START_REF] Huck | Minimum dispersion constrained nonnegative matrix factorization to unmix hyperspectral data[END_REF], [START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[END_REF]. Its underlying premise is that each pixel/spectrum is a linear combination of the endmembers. To be physically interpretable, two constraints are imposed: the abundance nonnegativity constraint (ANC) and the abundance sum-to-one constraint (ASC) for each pixel [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. Considering both constraints, the fullyconstrained least-squares method (FCLS) was developed in [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF]. Another well-known algorithm is SUnSAL (for Sparse Unmixing by variable Splitting and Augmented Lagrangian). It addresses the same optimization problem by taking advantage of the alternating direction method of multipliers (ADMM) [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF].

All these unmixing algorithms hugely suffer from noisy data and outliers within bands. Indeed, in real hyperspectral images for remote sensing, a considerable proportion (about 20%) of the spectral bands are noisy with low SNR, due to the atmospheric effect such as water absorption [START_REF] Zelinski | Denoising hyperspectral imagery and recovering junk bands using wavelets and sparse approximation[END_REF]. These bands need to be removed prior to applying any existing unmixing method; otherwise, the unmixing quality drastically decreases. Such sensitivity to outliers is due to the investigated 2 -norm as a cost function in the FCLS and SUnSAL algorithms, as well as all unmixing algorithms that explore least-squares solutions.

Information theoretic learning provides an elegant alternative to the conventional minimization of the 2 -norm in least-squares problems, by considering the maximization of the so-called correntropy [START_REF] Liu | Correntropy: properties and applications in non-gaussian signal processing[END_REF], [START_REF] Principe | Information theoretic learning: Renyi's entropy and kernel perspectives[END_REF], [START_REF] Wu | Robust hammerstein adaptive filtering under maximum correntropy criterion[END_REF], [START_REF] Chen | Convergence of a fixed-point algorithm under maximum correntropy criterion[END_REF]. Due to its stability and robustness to noise and outliers, the correntropy maximization is based on theoretical foundations and has been successfully applied to a wide class of applications, including cancer clustering [START_REF] Wang | Non-negative matrix factorization by maximizing correntropy for cancer clustering[END_REF], face recognition [START_REF] He | Maximum correntropy criterion for robust face recognition[END_REF], and recently hyperspectral unmixing [START_REF] Wang | Robust hyperspectral unmixing with correntropy-based metric[END_REF], to name a few. In these works, the resulting problem is optimized by the half-quadratic technique [START_REF] Nikolova | Analysis of half-quadratic minimization methods for signal and image recovery[END_REF], either in a supervised manner [START_REF] He | Maximum correntropy criterion for robust face recognition[END_REF] or as an unsupervised nonnegative matrix factorization [START_REF] Wang | Non-negative matrix factorization by maximizing correntropy for cancer clustering[END_REF], [START_REF] Wang | Robust hyperspectral unmixing with correntropy-based metric[END_REF].

In this paper, we derive the the fully-constrained (ANC and ASC) correntropy maximization problem for unmixing hyperspectral images. We propose to solve this problem with the ADMM algorithm. Indeed, the ADMM approach splits a hard problem into a sequence of small and handful ones [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. Its relevance to solve nonconvex problems was studied in [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]Section 9]. We show that ADMM provides a relevant framework for incorporating different constraints raised in unmixing problem. The resulting algorithm is called CUSAL for Correntropy-based Unmixing by variable Splitting and Augmented Lagrangian. The rest of the paper is organized as follows. We first revisit the classical unmixing problems with the linear model in Section II. In Section III, we present the proposed correntropybased unmixing problem. The resulting optimization problem is solved by the ADMM in Section IV. Experiments on synthetic and real images are presented in Section V and VI, respectively. Finally, Section VII provides some conclusions and future works.

II. CLASSICAL UNMIXING PROBLEM

The linear mixture model (LMM) assumes that each spectrum can be expressed as a linear combination of a set of pure material spectra, termed endmembers [START_REF] Keshava | Spectral unmixing[END_REF]. Consider a hyperspectral image and let Y ∈ R L×T denote the matrix of the T pixels/spectra of L spectral bands. Let y * t be its t-th column and y l * its l-th row, representing the l-th band of all pixels. For notation simplicity, we denote y t = y * t , for t = 1, 2, . . . , T . The LMM can be written as

y t = R r=1 x rt m r + n t = M x t + n t , (1)
where

M = [m 1 m 2 • • • m R ] ∈ R L×R is the matrix composed by the R endmembers, x t = [x 1t x 2t • • • x Rt ]
is the abundance vector associated with the t-th pixel, and n t ∈ R L is the additive noise. In matrix form for all pixels, we have

Y = M X + N , where X = [x 1 x 2 • • • x T ] ∈ R R×T
and N is the additive noise matrix.

In the following, the endmembers are assumed known, either from ground-truth information or by using any endmember extraction technique [START_REF] Honeine | Geometric unmixing of large hyperspectral images: a barycentric coordinate approach[END_REF]. The spectral unmixing problem consists in estimating the abundances for each pixel, often by solving the least-squares optimization problem

min xt y t -M x t 2 2 ,
for each t = 1, 2, . . . , T , where • 2 denotes the conventional 2 -norm. To be physically interpretable as emphasized in [START_REF] Keshava | Spectral unmixing[END_REF], the fractional abundances should be nonnegative (ANC) and satisfy the sum-to-one constraint (ASC). Considering both constraints, the fully-constrained least-squares problem [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF] is formulated as

min xt y t -M x t 2 2 , subject to x t ≥ 0 and 1 x t = 1, for t = 1, . . . , T,
where 1 ∈ R R×1 denotes the column vector of ones. Since there is no closed-form solution when dealing with the nonnegativity constraint, several iterative techniques have been proposed, such as the active set scheme with the Lawson and Hanson's algorithm [START_REF] Lawson | Solving Least Squares Problems[END_REF], the multiplicative iterative strategies [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF], and the fully-constrained least-squares (FCLS) technique [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF]. More recently, the alternating direction method of multipliers (ADMM) has been applied with success for hyperspectral unmixing problem, with the so-called SUnSAL algorithm [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF].

On sensitivity to outliers

All the aforementioned algorithms rely on solving a (constrained) least-squares optimization problem, thus inheriting the drawbacks of using the 2 -norm as the fitness measure. A major drawback is its sensitivity to outliers, where outliers are some spectral bands that largely deviate from the rest of the bands. Indeed, considering all the image pixels, the leastsquares optimization problems take the form

min X L l=1 y l * -(M X) l * 2 2 , ( 2
)
subject to any of the aforementioned constraints.

In this expression, (•) l * denotes the l-th row of its argument. From this formulation, it is easy to see how the squared 2 -norm gives more weight to large residuals, namely to outliers in which predicted values (M X) l * are far from actual observations y l * . Moreover, it is common for hyperspectral images to present up to 20% of unusable spectral bands due to low signal-to-noise ratio essentially from atmospheric effects, such as water absorption. In the following section, we propose to overcome this difficulty by considering the correntropy maximization principle from the information theoretic learning, which yields an optimization problem that is robust to outliers.

III. THE CORRENTROPY-BASED UNMIXING PROBLEM

In this section, we examine the correntropy and rewrite the hyperspectral unmixing problem as a correntropy maximization one.

A. Correntropy

The correntropy, proposed in [START_REF] Liu | Correntropy: properties and applications in non-gaussian signal processing[END_REF], [START_REF] Principe | Information theoretic learning: Renyi's entropy and kernel perspectives[END_REF], is a nonlinear, local similarity measure. For two random variables U and V , the correntropy is defined by

IE[κ(U, V )], (3)
where IE[•] is the expectation operator, and κ(•, •) is a shift-invariant kernel satisfying the Mercer theorem [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. In practice, while the joint distribution function of U and V is often unavailable, the sample estimator of correntropy is adopted instead.

Employing a finite number of data {(u l , v l )} L l=1 , the latter is estimated by

1 L L l=1 κ(u l , v l ), (4)
up to a normalization factor. Regarding the kernel function, the Gaussian kernel is the most commonly-used kernel for correntropy [START_REF] Liu | Correntropy: properties and applications in non-gaussian signal processing[END_REF], [START_REF] He | Maximum correntropy criterion for robust face recognition[END_REF], [START_REF] Chen | Maximum correntropy estimation is a smoothed map estimation[END_REF]. It is defined by

κ(u, v) = exp( -1 2σ 2 u -v 2 )
, where σ denotes the kernel bandwidth. This leads to the following expression for the correntropy:

1 L L l=1 exp -1 2σ 2 u l -v l 2 2 . ( 5
)
The maximization of the similarity measured by the correntropy, given by max 1 L L l=1 κ(u l , v l ), is termed the maximum correntropy criterion (MCC) [START_REF] Liu | Correntropy: properties and applications in non-gaussian signal processing[END_REF], [START_REF] Chen | Steady-state mean-square error analysis for adaptive filtering under the maximum correntropy criterion[END_REF]. It is noteworthy that well-known second-order statistics, such as the mean square error (MSE) depends heavily on the Gaussian and linear assumptions [START_REF] Liu | Correntropy: properties and applications in non-gaussian signal processing[END_REF]. However, in presence of non-Gaussian noise and in particular large outliers, i.e., observations greatly deviated from the data bulk, the effectiveness of the MSE-based algorithms will significantly deteriorate [START_REF] Wu | Robust hammerstein adaptive filtering under maximum correntropy criterion[END_REF]. By contrast, the MCCbased algorithms is appropriate for non-Gaussian signal processing, and is robust in particular against large outliers, as shown next.

We examine the behavior of the correntropy in terms of the residual error defined by l = u lv l . Therefore, the correntropy (5) becomes

1 L L l=1 exp -1 2σ 2 2 l .
Compared with the secondorder statistics, e.g. MSE, the correntropy is more robust with respect to the outliers, as shown in Fig. 1 illustrating the second-order and the correntropy objective functions in terms of the residual error. As the residual error augments, the secondorder function keeps increasing dramatically. On the contrary, the correntropy is only sensitive within a region of small residual errors, this region being controlled by the kernel bandwidth. For large magnitudes of residual error out of the central range, the correntropy grows extremely gently and tends to a constant. Consequently, the correntropy criterion is robust to large outliers, especially for small kernel bandwidths.

B. Correntropy-based unmixing problem

The correntropy-based unmixing problem consists in estimating the unknown abundance matrix X, by minimizing the objective function C (the negative of correntropy), given by

C(X) = - L l=1 exp -1 2σ 2 y l -(M X) l * 2 2 , ( 6
)
where the Gaussian kernel is considered and the nonnegative constraint is enforced on all the entries of X. Considering both the ANC and ASC constraints, the fully-constrained case for problem [START_REF] Nguyen | Supervised nonlinear unmixing of hyperspectral images using a preimage method[END_REF] becomes

min X C(X),
subject to X ≥ 0 and 1 x t = 1, for t = 1, . . . , T.

IV. ADMM FOR SOLVING THE CORRENTROPY-BASED UNMIXING PROBLEM

We first briefly review the alternating direction method of multipliers (ADMM), following the expressions in [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]Chap. 3]. Assuming that the functions f and g are closed, proper and convex, the ADMM solves the problem min f (x) + g(z) subject to Ax + Bz = c, [START_REF]Estimating abundance fractions of materials in hyperspectral images by fitting a post-nonlinear mixing model[END_REF] with variables x ∈ R n and z ∈ R m , where A ∈ R p×n and B ∈ R p×m and c ∈ R p . The Algorithm 1 The Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] Input: functions f and g, matrices A and B, vector c, parameter ρ 1: Initialize x 0 , z 0 and u 0 2: repeat 3:

x k+1 = arg min x f (x) + ρ 2 Ax + Bz k -c + u k 2 2 ; 4: z k+1 = arg min z g(z)+ ρ 2 Ax k+1 +Bz -c+ u k 2 2 ;
5:

u k+1 = u k + Ax k+1 + Bz k+1 -c; 6:
until stopping criterion algorithm ADMM, in its scaled form, is summarized in Algorithm 1.

In the following, we apply the ADMM algorithm to solve the correntropy-based unmixing problem in the fully-constrained case, presented in [START_REF] Chen | Nonlinear estimation of material abundances of hyperspectral images with 1 -norm spatial regularization[END_REF]. The main steps are summarized in Algorithm 2. Rewrite the variables to be optimized in a vector x ∈ R RT ×1 , which is stacked by the columns of the matrix

X, namely x = [x 1 x 2 • • • x T ] . Rewrite the vectors z = [z 1 z 2 • • • z T ] ∈ R RT ×1 and u = [u 1 u 2 • • • u T ] ∈ R RT ×1 , where z t = [z 1t z 2t • • • z Rt ] and u t = [u 1t u 2t • • • u Rt ] ,
for t = 1, 2, . . . , T . Following the formulation in Algorithm 1, we set

f (x) = C(x) + T t=1 ι {1} (1 x t ) (9) g(x) = ι R RT + (x) (10)
A = -I, B = I and c = 0, [START_REF] Huck | Minimum dispersion constrained nonnegative matrix factorization to unmix hyperspectral data[END_REF] where I is the identity matrix, 0 ∈ R RT ×1 is the zero vector and ι S (x) denotes the indicator function of the set S defined by

ι S (x) = 0, if x ∈ S ∞, elsewise.
In this case, the subproblem of the x-update in line 3 of Algorithm 1 addresses a nonconvex problem without closed-form solution. To overcome this difficulty, we apply an inexact ADMM variant in lines 3-5 of Algorithm 2, which solves the subproblem iteratively using the gradient descent method, instead of solving it exactly and explicitly. Before that, we eliminate the T equality constraints, i.e., the sum-to-one constraints, by replacing x Rt with x Rt = 1 -R-1 r=1 x rt , for t = 1, 2, ..., T . Let x ∈ R (R-1)T ×1 be the reduced vector of variables, stacked by

x t = [x 1t x 2t • • • x (R-1)t ] , t = 1,
2, ..., T . By this means, the function ( 9) is transformed into the reduced-form

f 1 (x) = L l=1 -exp -T t=1 h l (x t ) 2 2σ 2 , ( 12
)
where

h l (x t ) = y lt -m lR -R-1 p=1 (m lp -m lR )
x pt , for l = 1, 2, ..., L. The gradient of ( 12) with respect to x is stacked as

∂f 1 ∂x = [ ∂f 1 ∂x 1 • • • ∂f 1 ∂x T ] ∈ R (R-1)T ×1 , where ∂f 1 ∂xt = [ ∂f 1 ∂x 1t • • • ∂f 1 ∂x (R-1)t ]
, with the entries given by

∂f 1 (x) ∂x rt = 1 σ 2 L l=1 (m lR -m lr ) exp -T s=1 h l (xs) 2 2σ 2 h l (x t ),
for all r = 1, 2, ..., (R -1) and t = 1, 2, ..., T . Similarly, the function ρ

2 x-z k -u k 2 2 is expressed with respect to x as φ(x)= ρ 2 T t=1 R-1 r=1 (x rt -z rt,k -u rt,k ) 2 + 1- R-1 r=1 x rt -z Rt,k -u Rt,k 2 ,
with the entries in its gradient ∂φ ∂x given by

∂φ(x) ∂x rt = ρ(x rt + R-1 p=1 x pt -1+z Rt,k -z rt,k +u Rt,k -u rt,k ), (13)
for all r = 1, 2, ..., (R -1) and t = 1, 2, ..., T .

The solution of the z-update in line 4 Algorithm 1 is the projection of x k+1 -u k onto the first orthant as shown in line 7 of Algorithm 2.

V. EXPERIMENTS WITH SYNTHETIC DATA

In this section, the performance of the proposed CUSAL is evaluated on synthetic data. A comparative study is performed considering five state-of-the-art methods proposed for both linear and nonlinear unmixing models. We apply a three-fold stopping criterion for Algorithm 2 according to [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF]: (i) the primal and dual residuals are small enough, (ii) the primal residual starts to increase, or (iii) the predefined maximum Algorithm 2 Correntropy unmixing (CUSAL)

1: Initialize k = 0, ρ > 0, η > 0, σ > 0, norms of the primal and dual residuals r 0 = 10 5 and s 0 = 10 5 ; x 0 , z 0 and u 0 ; 2: repeat 3: repeat 4:

x k+1 = x k+1 -η( ∂f 1 ∂x k+1 + ∂φ ∂x k+1 ); 5:
until convergence 6: reform x k+1 using x k+1 ;

7:

z k+1 = max(0, x k+1 -u k ); 8: u k+1 = u k -(x k+1 -z k+1 ); 9: r k+1 = x k+1 -z k+1 2 , s k+1 = ρ z k+1 - z k 2 ;
10:

k = k + 1; 11: until stopping criterion iteration number is attained.
The proposed CUSAL algorithm is compared to the following methods:

• Fully-Constrained Least-Squares (FCLS) [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF]:

The FCLS is developed for the linear model. Enforcing both ANC and ASC constraints, this technique yields the optimal abundance matrix in the least-squares sense. See also [START_REF] Honeine | Geometric unmixing of large hyperspectral images: a barycentric coordinate approach[END_REF]. • Sparse Unmixing by variable Splitting and Augmented Lagrangian (SUnSAL) [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF]: This method is based on the ADMM. Under the linear assumption, several variants are developed by including different constraints. In the following, we consider the fully-constrained variant, termed SUnSAL-FCLS. • The Bayesian algorithm for Generalized Bilinear Model (BayGBM) [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], [START_REF]Unmixing hyperspectral images using the generalized bilinear model[END_REF]: The BayGBM is a nonlinear abundance estimation method based on the generalized bilinear model (GBM), which adds the second-order interactions between endmembers to the linear model and is defined as

y t = M x t + R-1 i=1 R j=i+1 γ ij,t x it x jt (m i m j )+n t ,
where 0 ≤ γ ij,t ≤ 1 controls the interactions between endmembers m i and m j , and is the Hadamard product (element-wise product) operator. The BayGBM is designed considering both ANC and ASC.

• The Bayesian algorithm for Polynomial Postnonlinear Mixing Model (BayPPNMM) [START_REF] Altmann | Supervised nonlinear spectral unmixing using a postnonlinear mixing model for hyperspectral imagery[END_REF]:

The BayPPNMM is a Bayesian-based algorithm proposed to estimate the parameters involved in the polynomial postnonlinear mixing model (PPNMM). This model assumes that the pixel reflectances are nonlinear functions of endmembers, and is defined by

y t = M x t + b t (M x t ) (M x t ) + n t , ( 14
)
where the nonlinear terms are characterized by b t ∈ R, and both ANC and ASC are required.

• Kernel

Fully-Constrained Least-Squares (KFCLS) [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF]: This method generalizes FCLS, by replacing the inner product with a kernel function. In the following, the Gaussian kernel is applied for simulation. A preliminary test on 100 randomly selected pixels is performed for each data under study, using the set {1, 2, 3, 4, 5} for the parameter σ and the set {1, 10 -1 , 10 -2 , 10 -3 , 10 -4 } for the parameter that controls the relaxation of the sum-to-one constraint.

A series of experiments are performed, mainly considering the influences of three aspects: (i) noise level, (ii) number of corrupted bands and (iii) number of endmembers.

In the following simulations, each synthetic image has a fixed size of 50 × 50 pixels. First, it is generated using the linear mixture model [START_REF] Keshava | Spectral unmixing[END_REF], where Gaussian noise n t of SNR ∈ {15, 35}dB is added. The R ∈ {3, 6, 9} endmembers are randomly drawn from 20 spectra selected from the United States Geological Survey (USGS) digital spectral library [START_REF] Bioucas-Dias | Hyperspectral subspace identification[END_REF]. These candidate endmembers are measured over L = 244 continuous bands with the wavelength ranging from 0.2μm to 3.0μm. The abundance vectors x t are uniformly generated using a Dirichlet distribution as in [START_REF] Bioucas-Dias | Hyperspectral subspace identification[END_REF], [START_REF] Halimi | Unsupervised unmixing of hyperspectral images accounting for endmember variability[END_REF]. To imitate the noisy bands in the real hyperspectral images, several bands in the generated data are corrupted by replacing the corresponding rows y l * with random values within [0, 1]. The number of corrupted bands varies in the set {0, 20, 40, 60}.

The unmixing performance is evaluated using the abundance root mean square error (RMSE) [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], [START_REF] Yokoya | Nonlinear unmixing of hyperspectral data using semi-nonnegative matrix factorization[END_REF], defined by

RMSE = 1 RT T t=1 x t -x t 2 ,
where x t is the estimated abundance vector. Fig. 2 illustrates the average of RMSE over 10 Monte-Carlo realizations. It is easy to see that, in presence of outlier bands, the proposed CUSAL outperforms all the comparing methods in terms of RMSE, for different noise levels and numbers of endmembers.

It is also shown that the performance of the proposed algorithm improves when increasing the SNR.

VI. EXPERIMENTS WITH REAL DATA

This section presents the performance of the proposed algorithm on a real hyperspectral image. We consider a 250 × 190 sub-image taken from the Cuprite mining image, acquired by the AVIRIS sensor when flying over Las Vegas, Nevada, USA. The image has been widely investigated in the literature [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF], [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF]. The raw data contains L = 224 bands, covering a wavelength range 0.4 -2.5μm. Among, there are 37 relatively noisy ones with low SNR, namely the bands 1 -3, 105 -115, 150 -170, and 223 -224. The geographic composition of this area is estimated to include 14 minerals [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF]. Neglecting the similar signatures, a reduced number of 12 endmembers is considered [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF], [START_REF] Lu | Manifold regularized sparse nmf for hyperspectral unmixing[END_REF] The unmixing results in terms of the averaged spectral angle distance (SAD) between the input spectra y t and the reconstructed ones y t are illustrated in Fig. 3. The SAD is defined by

SAD = 1 T T t=1
arccos y t y t y t y t .

In absence of noisy bands (i.e., L = 187 bands), all the considered methods lead to satisfactory abun- dance maps, with BayPPNMM providing the smallest SAD. As the number of noisy bands increases, especially from L = 199 to L = 224, the unmixing performance of the state-of-the-art methods deteriorates drastically, while the proposed CUSAL yields stable SAD. The obtained results confirm the good behavior of the proposed CUSAL and its robustness in presence of corrupted spectral bands.

VII. CONCLUSION

This paper presented a supervised unmixing algorithm based on the correntropy maximization principle. The alternating direction method of multipliers (ADMM) was investigated in order to maximize the correntropy criterion under nonnegativity and sum-to-one contraints. The effectiveness and robustness of the proposed method were validated on both synthetic and real hyperspectral images. Future works include the generalization of the correntropy criterion to account for the multiple reflection phenomenon, as well as a sparsity-aware version of the proposed method.

Fig. 1 .

 1 Fig. 1. Illustration of second-order objective function (solid line) and the correntropy objective function (dashed lines) with different values of the kernel bandwidth.

  . The VCA technique is first applied to extract the endmembers on the clean image with L = 187 bands. The resulting 12 pixel indices are retained. Starting from L = 187 bands, the noisy bands, randomly chosen from the bands 1 -3, 105 -115, 150 -170, and 223 -224, are gradually included to form a series of input data. In practice, the experiments are conducted with L = 187, 193, 199, 205, 211, 217, 223 and 224 bands.

9 Fig. 2 .

 92 Fig. 2. The averaged root mean square error (RMSE) with respect to the number of corrupted bands (10 Monte-Carlo realizations). Left to right: Various number of endmembers. Top to bottom: Different SNR.

Fig. 3 .

 3 Fig. 3. The SAD values using different numbers of bands on the Cuprite image.

ACKNOWLEDGMENT

This work was supported by the French ANR, grant HYPANEMA: ANR-12BS03-0033.


Abderrahim Halimi

Gerald Buller

Steve Mclaughlin

Paul Honeine

FILTERING SMOOTH ALTIMETRIC SIGNALS USING A BAYESIAN ALGORITHM

Keywords: Altimetry, Bayesian algorithm, coordinate descent algorithm, gamma Markov random fields

This paper presents a new Bayesian strategy for the estimation of smooth signals corrupted by Gaussian noise. The method assumes a smooth evolution of a succession of continuous signals that can have a numerical or an analytical expression with respect to some parameters. The Bayesian model proposed takes into account the Gaussian properties of the noise and the smooth evolution of the successive signals. In addition, a gamma Markov random field prior is assigned to the signal energies and to the noise variances to account for their known properties. The resulting posterior distribution is maximized using a fast coordinate descent algorithm whose parameters are updated by analytical expressions. The proposed algorithm is tested on satellite altimetric data demonstrating good denoising results on both synthetic and real signals. The proposed algorithm is also shown to improve the quality of the altimetric parameters when combined with a parameter estimation strategy.

INTRODUCTION

A satellite altimeter is a nadir-viewing radar that emits regular pulses and records the travel time, the magnitude and the shape of each return signal after reflection on the Earth's surface. This reflected echo provides information about some physical parameters such as the range between the satellite and the observed scene (denoted by τ ), the significant wave height (denoted by SWH) and the wind speed (related to the signal's amplitude P u ). The oceanic altimetric signals can be described as a succession of continuous functions corrupted by noise [START_REF] Hapke | Bidirectional reflectance spectroscopy. I. Theory[END_REF][START_REF] Brown | The average impulse response of a rough surface and its applications[END_REF][START_REF] Halimi | Parameter estimation for peaky altimetric waveforms[END_REF] (see Fig. 1 (top) that shows a succession of 800 signals acquired by the Jason-2 mission). Many recent studies and missions have focused on improving the quality of these signals by reducing the noise effect. This goal is achieved by amending the estimation algorithm to incorporate the known smooth properties of the altimetric parameters [START_REF] Maus | Improved ocean-geoid resolution from retracked ERS-1 satellite altimeter waveforms[END_REF][START_REF] Sandwell | Retracking ERS-1 altimeter waveforms for optimal gravity field recovery[END_REF][START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF] or by operating on the observed signals to reduce their noise [START_REF] Ollivier | Nouvelle approche pour l'extraction de paramètres géophysiques à partir des mesures en altimétrie radar[END_REF][START_REF] Thibaut | Singular value decomposition applied on altimeter waveforms[END_REF]. This paper adopts this second approach by exploiting the smooth evolution of the successive altimetric signals. The main motivation for this choice is to propose a denoising algorithm that is independent from the parameter estimation algorithm, thus, it can be easily combined with any existing estimation algorithms [START_REF] Maus | Improved ocean-geoid resolution from retracked ERS-1 satellite altimeter waveforms[END_REF][START_REF] Sandwell | Retracking ERS-1 altimeter waveforms for optimal gravity field recovery[END_REF][START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF][START_REF] Amarouche | Improving the Jason-1 ground retracking to better account for attitude effects[END_REF] leading to an improvement in parameter estimation. The first contribution of this paper is the development of a hierarchical Bayesian model to denoise a set of smooth signals. Each signal is assumed to be corrupted by additive, independent and non identically Gaussian noise. This noise model generalizes the independent and identically distributed (i.i.d.) Gaussian noise that is generally assumed when considering altimetric data [START_REF] Amarouche | Improving the Jason-1 ground retracking to better account for attitude effects[END_REF][START_REF] Halimi | A semianalytical model for delay/Doppler altimetry and its estimation algorithm[END_REF]. A gamma Markov random field (GMRF) prior [START_REF] Dikmen | Gamma markov random fields for audio source modeling[END_REF] is considered to account for the correlation between the noise variances and to better approximate the speckle noise. The signal energies are also assigned a GMRF prior to better approximate their continuity. Note that the proposed Bayesian hierarchy is generic in the sense that it does not assume a specific signal model. Indeed, the signal can be expressed by a numerical formula or by a given linear/nonlinear analytical function with respect to (w.r.t.) some parameters.

The second contribution of this paper is the derivation of a denoising algorithm associated with the proposed hierarchical Bayesian model. The minimum mean square error (MMSE) and maximum a posteriori (MAP) estimators of the unknown signals/parameters cannot be easily computed from the obtained joint posterior. In this paper, the MAP estimator is evaluated by considering a coordinate descent algorithm (CDA) [START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF][START_REF] Bertsekas | Nonlinear programming[END_REF]. The resulting algorithm sequentially updates the noiseless signals, the noise variances and other hyperparameters by analytical formulas leading to a reduced computational cost. The proposed Bayesian model and estimation algorithm are validated using synthetic and real altimetric data acquired during the Jason-2 mission. The obtained results are very promising and show the potential of the proposed denoising strategy.

The paper is structured as follows. Section II introduces the observation model and the considered altimetric signal. The proposed hierarchical Bayesian model and its estimation algorithm are introduced in Sections III and IV. Section V validates the proposed technique using simulated data with controlled ground truth. Section VI shows results obtained using real data obtained from the Jason-2 mission. Finally, conclusions and future work are reported in Section VII.

PROBLEM FORMULATION

Consider M successive signals S ∈ R K×M and let Y ∈ R K×M denote their noisy version. Let y :m ∈ R K×1 be the mth column of Y and y k: ∈ R 1×M its kth row, representing the kth temporal gate for all signals. For notation simplicity, we denote y :m = y m , ∀m, and y k: = y k , ∀k. Using these notations, the observation model is given by

y m = s m (Θ m ) + e m , with e m ∼ N (0 K , Σ) (1)
where y m and s m are (K × 1) vectors representing the mth observed and noiseless signals, 0 K denotes the (K × 1) vector of zeros and e m is a centered Gaussian noise vector with a diagonal covariance matrix

Σ = diag σ 2 with σ 2 = σ 2 1 , • • • , σ 2 K T a (K × 1) vector. The signals S
might depend on some parameters (by a linear or nonlinear expression) which are denoted by the (1 × H) vector

Θ m = [θ 1 (m), • • • , θ H (m)]
containing the H parameters of the mth signal. Note, however, that the proposed method does not necessarily require a parametric expression for S, and is valid provided that the signals satisfy some properties (as described in the following). The smooth variation of the oceanic altimetric signals (see Fig. 1 (bottom)) can be highlighted by expressing the model (1) as follows

y k = s k (Θ) + e k , with e k ∼ N 0 M , σ 2 k I M (2)
where k ∈ {1, • • • , K} indexes the signal samples that are known as "temporal gates", I M denotes the (M × M ) identity matrix and s k is a smooth (M × 1) vector representing the signal evolution at the kth gate (see Fig. 1 (bottom) for examples). The proposed Bayesian method aims to filter the observed signals

y k , k ∈ {1, • • • , K} to retrieve the noiseless signals s k , k ∈ {1, • • • , K}.

Conventional altimetric model

This section introduces the satellite altimetric model that will be considered in this paper since it satisfies the model described above. The altimetric nonlinear model accounts for three parameters (P u , τ, SWH) as follows [START_REF] Brown | The average impulse response of a rough surface and its applications[END_REF][START_REF] Amarouche | Improving the Jason-1 ground retracking to better account for attitude effects[END_REF]

s(t) = P u 2 1 + erf t -τ s -ασ 2 c √ 2σ c exp -α t -τ s - ασ 2 c 2 (3) where σ 2 c = SWH 2c 2 + σ 2 p , erf (t) = 2 √ π t 0 e -z 2 dz is the Gaussian error function, t is time, τ s = 2τ
c (resp. τ ) is the epoch expressed in seconds (resp. meters), c is the speed of light, α and σ 2 p are two known parameters (depending on the satellite and on the measurement instrument). Note that the discrete altimetric signal is gathered in the vector s

= (s 1 , • • • , s K ) T , where K = 104 gates, s k = s (kT ), T is the time resolution and Θ m = [SWH(m), τ (m), P u (m)] is a (1× 3
) vector containing the 3 altimetric parameters SWH, τ, P u for the mth signal.

The altimetric signals are corrupted by speckle noise that, thanks to the averaging that takes place on-board the satellite, can be approximated by additive Gaussian noise as shown in [START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF][START_REF] Martin-Puig | SAR altimeter retracker performance bound over water surfaces[END_REF][START_REF] Germain | A revisit to the GNSS-R code range precision[END_REF][START_REF] Halimi | Cramér-Rao bounds and estimation algorithms for delay/Doppler and conventional altimetry[END_REF]. Thus, the observation altimetric model satisfies [START_REF] Hapke | Bidirectional reflectance spectroscopy. I. Theory[END_REF]. Moreover, the noise variances obtained, σ 2 k , k ∈ {1, • • • , K}, after the satellite averaging, are correlated due to the nature of the speckle noise (this correlation will be considered in the proposed Bayesian scheme). Finally, it is worth-noting that this paper only considers oceanic observations which generally show a smooth variation between successive signals.

HIERARCHICAL BAYESIAN MODEL

This section introduces a hierarchical Bayesian model to estimate the unknown (K ×M ) matrix of noiseless signals S and the (K ×1) vector σ containing the noise variances associated with the M signals considered.

Likelihood

The observation model defined in (2) and the Gaussian properties of the noise sequence e

k , k ∈ {1, • • • , K}, yield f (y k |s k , σ 2 k ) ∝ 1 σ 2 k M 2 exp - ||y k -s k || 2 2σ 2 k ( 4
)
where ∝ means "proportional to", ||•|| denotes the standard l 2 norm such that ||x|| 2 = x T x and s k (Θ) has been denoted by s k for brevity. Assuming independence between the observed signals leads to f (Y |Θ, S)

∝ K k=1 f (y k |s k , σ 2 k ).

Priors for the observed signal

As previously assumed, the successive observed signals evolve slowly leading to smooth vectors s k , for k ∈ 1, • • • , K (see Fig. 1 (bottom)). This property is satisfied by considering a Gaussian prior for s k ensuring smoothness as follows

s k | 2 k ∼ N 0 M , 2 k H , (5)
where H is an (M × M ) matrix representing the squaredexponential covariance function given by H(m, m

) = exp - (m-m ) 2 (30) 2
, which introduces the correlation between the successive signals and 2 k is a variance parameter that is gate dependent. From [START_REF] Sandwell | Retracking ERS-1 altimeter waveforms for optimal gravity field recovery[END_REF], it is clear that this variance is related to the energy of the signals at the kth gate (via the norm s T k H -1 s k ). Moreover, because of the continuity of the signal s m w.r.t. the temporal gates, the signal energies vary smoothly from one gate to another. Therefore, we expect 2 k to vary smoothly from one gate to another which will be introduced by considering a specific prior for 2 k , as explained in Section III-D.

Prior for the noise variance

Due to the speckle origins of the corruption noise, we expect the noise variances σ 2 k , k ∈ {1, • • • , K} to vary smoothly. This behavior is considered by introducing an auxiliary vector w (of size K × 1) and assigning a gamma Markov random field prior (GMRF) for (σ, w) (see [START_REF] Dikmen | Gamma markov random fields for audio source modeling[END_REF][START_REF] Halimi | Bayesian filtering of smooth signals: Application to altimetry[END_REF] for more details regarding this prior). This prior ensures that each σ 2 k is connected to two neighboring elements of w and vice-versa via a coupling parameter ζ > 1 that controls the amount of correlation enforced by the GMRF (see Fig. 2 (a)). Note that the variances σ 2 k , σ 2 k for k = k are conditionally independent and that the correlation is introduced via the auxiliary variables w. An interesting property of this joint prior is that the conditional prior distributions of σ and w reduce to conjugate inverse gamma (IG) and gamma (G) distributions, respectively, as follows [START_REF] Dikmen | Gamma markov random fields for audio source modeling[END_REF]

σ 2 k |w k-1 , w k , ζ ∼ IG [2ζ, ζ (w k-1 + w k )] σ 2 K |w K-1 , ζ ∼ IG (ζ, ζw K-1 ) w 2 k |σ 2 k , σ 2 k+1 , ζ ∼ G 2ζ, ζ σ 2 k + ζ σ 2 k+1 -1 (6)
where k ∈ {1, • • • , K -1}.

Hyperparameter priors

As previously explained, the hyperparameters 2 k are closely related to the energies s T k H -1 s k . Considering this property and the continuity of the signal suggest the presence of a correlation between the parameters 2 k . This correlation can be introduced using a similar GMRF prior as in Section III-C for the couple ( , v), where v is a new auxiliary vector and η is the coupling parameter (see Fig. 2 (b)).

σ 2 2 w 2 • • • σ 2 K w 1 σ 2 1 w 0 ζ ζ ζ ζ ζ ζ (a) GMRF for σ 2 2 v 2 • • • 2 K v 1 2 1 v 0 η η η η η η (b) GMRF for

Posterior distributions

The parameters of interest are X = (S, σ, w, , v). The joint posterior distribution of this Bayesian model can be computed using the following hierarchical structure

f (X|Y ) ∝ f (Y |S, σ)f (S| ) f ( , v) f (σ, w) , (7)
where we have assumed a priori independence between the parameters. For simplicity, f (x|θ) has been denoted by f (x) when the parameter θ is a user-fixed parameter.

COORDINATE DESCENT ALGORITHM

This section describes the optimization algorithm maximizing the posterior (7) w.r.t. X which provides the parameter MAP estimator. Because of the large number of parameters in X = (S, σ, w, , v), we propose a coordinate descent algorithm [START_REF] Bertsekas | Nonlinear programming[END_REF] that sequentially updates the different parameters. More precisely, in each step, the posterior distribution is maximized w.r.t. one parameter, the others being fixed. This process is repeated until the algorithm has converged to a local minimum of the negative log-posterior C(S, σ, w, , v). Thus, the algorithm iteratively updates each parameter by maximizing its conditional distribution (similarly to the iterated conditional modes algorithm [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF]) which are given by

• s k |y k , σ 2 k , 2 k ∼ N (s k , Γ k ) • σ 2 k |y k , s k , w k ∼ IG 2ζ + M 2 , β1 2 • 2 k |y k , s k , v k ∼ IG 2η + M 2 , β2 2 • w 2 k |σ 2 k , σ 2 k+1 , ζ ∼ G 2ζ, ζ σ 2 k + ζ σ 2 k+1 -1 • v 2 k | 2 k , 2 k+1 , η ∼ G 2η, η 2 k + η 2 k+1 -1
where

s k = 1 σ 2 k Γ k y k ,, Γ k = H -1 2 k + I M σ 2 k -1 , β 1 = ||y k - s k || 2 +2ζ(w k-1 +w k ) and β 2 = s T k H -1 s k +2η(v k-1 +v k ).
The mode of each conditional distribution is uniquely and analytically obtained which reduces the computational cost of the algorithm. In this paper, the parameters have been initialized as follows:

σ (0) = s (0) m = 1 M M n=1 y n , ∀m, 2 k = 10, ∀k, and w (0) k = v (0) k = 10 -12
, ∀k. Note that more elaborate initialization procedures can be investigated, but these proposed values provided minima of good quality in the simulations considered (see Sections V and VI). More details regarding the stopping criteria of the algorithm are provided in [START_REF] Halimi | Bayesian filtering of smooth signals: Application to altimetry[END_REF].

SIMULATION RESULTS ON SYNTHETIC DATA

This section evaluates the performance of the proposed SSE (for smooth signal estimation) algorithm on synthetic data when considering two experiments. The first experiment evaluates the performance when varying the number of the denoised signals M . Indeed, the SSE algorithm considers successive sets of length M to denoise the observed altimetric signals. Therefore, N = 5000 signals are generated according to the altimetric model (3) while using a realistic variation of the altimetric parameters I. For example, for a length set M = 250, the algorithm is run 20 times to process the N = 5000 signals. The performance is evaluated using the reconstruction signal to noise ratio (RSNR) given by RSNR = 10 log 10

Θ m = [SWH(m), τ (m), P u (m)].
M m=1 ||sm|| 2 M m=1 ||sm-sm|| 2
, where s m (resp. s m ) denotes the noiseless (resp. denoised) signal. Overall, these results show an ≈ 11 dB improvement in the processed data with an increasing RSNR w.r.t. M . However, a high number of M requires higher computational cost (mainly due to the matrix inversion in (IV)), while too small M leads to more iterations. M = 500 represents a good compromise and we consider this value for the rest of the paper [START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF].

The second experiment compares our algorithm with the state-of-the-art singular value decomposition (SVD) filtering strategy proposed in [START_REF] Ollivier | Nouvelle approche pour l'extraction de paramètres géophysiques à partir des mesures en altimétrie radar[END_REF][START_REF] Thibaut | Singular value decomposition applied on altimeter waveforms[END_REF] (with a threshold equal to 84%). The study is performed when varying SWH ∈ [0.5, 8] m with fixed τ = 31 gates and P u = 130. For each SWH, 500 synthetic signals are generated using the model (3) with different noise realizations (500 Monte carlo runs) and processed using the two considered algorithms. Table II reports the obtained RSNR when considering SVD and SSE for different SWH. This table shows an average improvement by 6 dB when considering SVD and by 12 dB when considering the proposed SSE algorithm. These results highlight the potential of the proposed strategy in denoising the altimetric signals and improving the estimation of the altimetric parameters.

RESULTS ON JASON-2 REAL DATA

This section is devoted to the validation of the proposed SSE denoising algorithm when applied to the oceanic Jason-2 dataset. The data considered cover a period of 36 minutes and consist of 43000 real signals that were extracted from pass 30 of cycle 35. Figs. 1 (top) and (middle) present a sequence of 800 Jason-2 signals before and after filtering. These figure clearly show a reduction in the noise affecting the signals after the application of the SSE algorithm especially in the tail of the signal (the decreasing part), which was most affected by the speckle noise. Table III shows the effect of the denoising algorithms on the estimated altimetric parameters that are obtained using the well known least-squares (LS) estimation based strategy [START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF][START_REF] Amarouche | Improving the Jason-1 ground retracking to better account for attitude effects[END_REF]. This table compares the classical LS algorithm (applied to noisy data) with SVD-LS and the proposed SSE-LS (LS applied to filtered data). This table shows a good agreement between the means of the estimated parameters for the LS, SVD-LS and SSE-LS algorithms (except P u that is slightly reduced by SSE-LS). Moreover, the smallest parameter standard-deviations (STDs) are achieved by the proposed SSE-LS. This STD reduction is of great importance for many practical applications related to oceanography such as bathymetry. Comparing SSE-LS to LS, Table III highlights an STD improvement factor by 6 for SWH, 4 for τ and 5 for P u . Finally, Table III also compares the computational costs of the three considered algorithms when processing the 43000 signals (the result is reported for each signal). Because of the filtering step, both SVD-LS and SSE-LS require more computational time than LS. Note that the proposed SSE algorithm requires more computational time than the SVD approach. However, this cost (about 12% of additional computational times w.r.t. the LS algorithm) must be balanced by the performance improvement in terms of RSNR and parameter STDs. These results confirm the good performance of the proposed strategy for denoising smooth signals such as oceanic altimetric signals.

CONCLUSIONS

This paper has presented a new Bayesian strategy for denoising smooth signals corrupted by Gaussian noise. The method assumed a reduced evolution of successive continuous signals, which can have a numerical expression or be given by a linear/nonlinear function with respect to some parameters. A Bayesian model was proposed to take into account the Gaussian properties of the noise, and the smooth properties of the signal evolution. The resulting posterior distribution was maximized using a fast coordinate descent algorithm that showed good results on both synthetic and real altimetric signals. The proposed algorithm was also evaluated by combining it with a commonly used parameter estimation strategy for the altimetric parameters. The estimated parameters showed a clear improvement highlighting the benefit of the proposed algorithm. It is worth-noting that the proposed strategy is fast and generic and thus could be applied when considering other altimetric technologies such as delay/Doppler altimetry [START_REF] Halimi | A semianalytical model for delay/Doppler altimetry and its estimation algorithm[END_REF][START_REF] Raney | The delay/Doppler radar altimeter[END_REF][START_REF] Halimi | Including antenna mispointing in a semi-analytical model for delay/Doppler altimetry[END_REF]. This point will be considered in future work.

Fig. 1 .

 1 Fig. 1. (Top) Example of 800 noisy Jason-2 signals. (Middle) Filtered Jason-2 signals by the proposed SSE algorithm. (Bottom) Signal evolutions in gates (30, 50, 90) (continuous lines) and their smooth approximation with SSE (dashed lines).

Fig. 2 .

 2 Fig. 2. Proposed 1st order GMRF neighborhood structures for (a) the noise variances σ and (b) the signal energies .

  This realistic sequence of parameters is obtained by applying the CD-BM algorithm [6] on 5000 real Jason-2 signals (since it provides physically realistic smooth parameters) where we obtain SWH ∈ [3.4, 5.4] meters, τ ∈ [14.3, 15] meters and P u ∈ [150, 190] unit. The generated synthetic signals are then corrupted by speckle noise resulting from the averaging of L = 90 signals and leading to RSNR = 19.55 dB. The obtained N = 5000 signals are processed by the proposed algorithm while considering different set lengths as shown in Table

Table 1 .

 1 Performance of the proposed SSE algorithm w.r.t. the filter length (5000 signals). The corrupted data presents an RSNR = 19.55dB. Time per signal (ms) 0.35 0.24 0.25 0.38 1.03 5.00 17.47

			Filter length
		50	100 250 500 1000 2500 5000
	RSNR (dB)	31.1 31.4 31.5 31.6 31.7 31.7	31.7

Table 2 .

 2 RSNR (in dB) with respect to SWH. The corrupted data presents an RSNR = 19.55dB. SVD 26.35 26.43 26.30 26.02 26.03 26.07 26.08 25.92 25.86 SSE 32.24 32.21 32.22 32.13 32.15 32.10 32.22 32.13 32.07

					SWH (m)
	0.5	1	2	3	4	5	6	7	8
	RSNR (dB)

Table 3 .

 3 Performance on real Jason-2 data (45000 signals). SWH (cm) τ (cm) P u

		LS	242	14.68 167.73
	Mean	SVD-LS	241	14.67 166.62
		SSE-LS	248	14.68 164.83
		LS	59.9	12.01	6.18
	STD	SVD-LS	18.14	6.02	6.09
		SSE-LS	9.03	2.94	1.21
	Average time per signal (ms)	LS SVD-LS SSE-LS		8.56 9.05 9.63

This work was supported by the EPSRC Grants EP/JO15180/1, EP/N003446/1, and EP/K015338/1


Abderrahim Halimi

Paul Honeine

Jose Bioucas-Dias

Gerald S Buller

Steve Mclaughlin

NONLINEAR HYPERSPECTRAL UNMIXING ACCOUNTING FOR SPATIAL ILLUMINATION VARIABILITY

Keywords: Hyperspectral, nonlinear unmixing, Bayesian estimation, coordinate descent, gamma Markov random field

This paper presents a new supervised algorithm for nonlinear hyperspectral unmixing. Based on the residual component analysis model, the proposed model assumes the linear model to be corrupted by an additive term that accounts for bilinear interactions between the endmembers. The proposed formulation considers also the effect of the spatial illumination variability. The parameters of the proposed model are estimated using a Bayesian strategy. This approach introduces prior distributions on the parameters of interest to take into account their known constraints. The resulting posterior distribution is optimized using a coordinate descent algorithm which allows us to approximate the maximum a posteriori estimator of the unknown model parameters. The proposed model and estimation algorithm are validated on both synthetic and real images showing competitive results regarding the quality of the inferences and the computational complexity when compared to the state-of-the-art algorithms.

INTRODUCTION

Spectral unmixing (SU) consists of identifying the macroscopic materials (endmembers) present in a hyperspectral image and quantifying the proportions (abundances) of these materials in all pixels of the image. The linear mixture model (LMM) is the widely used model for SU mainly because of its simplicity. However, this model can be inappropriate for some hyperspectral scenarios, namely if there are volumetric scattering, or terrain relief, or intimate mixtures of materials [START_REF] Hapke | Bidirectional reflectance spectroscopy. I. Theory[END_REF]. Nonlinear (NL) mixture models (NLMM) appear then as an alternative to better account for those effects [START_REF] Heylen | A review of nonlinear hyperspectral unmixing methods[END_REF][START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF]. There exists two main families for NLMMs: the first family is signal processing based and seeks to construct flexible models that can represent a wide range of nonlinearities. The second family is physical based which includes models accounting for multiple scattering such as polynomial [START_REF] Altmann | Supervised nonlinear spectral unmixing using a postnonlinear mixing model for hyperspectral imagery[END_REF] or bilinear models [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF][START_REF]Unmixing hyperspectral images using the generalized bilinear model[END_REF][START_REF] Bioucas-Dias | Nonlinear mixture model for hyperspectral unmixing[END_REF][START_REF] Fan | Comparative study between a new nonlinear model and common linear model for analysing laboratory simulated-forest hyperspectral data[END_REF][START_REF] Meganem | Linear-quadratic mixing model for reflectances in urban environments[END_REF].

This paper introduces a residual component (RC) mixture model for nonlinear hyperspectral unmixing. The model generalizes the well-known LMM by accounting for the spatial illumination variation and the presence of a nonlinear residual term. This term is a modification of the polynomial term proposed in [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF], to account for the bilinear interactions between the endmembers. Estimating the abundances associated with this mixture model is a challenging problem. We propose here a hierarchical Bayesian model to estimate the parameters and hyperparameters associated with the RC model. This hierarchical model introduces prior distributions that enforce known physical constraints on the estimated parameters such as the sum-to-one and positivity of the abundances, and the positivity of the nonlinear coefficients. Moreover, the spatial correlation of the residual term has been introduced by considering Markov random fields (MRF) [START_REF] Dikmen | Gamma markov random fields for audio source modeling[END_REF]. Using the likelihood and the considered prior distributions, the joint posterior distribution of the unknown parameter vector is then derived. The minimum mean square error (MMSE) and maximum a posteriori (MAP) estimators of these parameters cannot be easily computed from the obtained joint posteriors. In this paper, the MAP estimator is evaluated by considering a coordinate descent algorithm (CDA) [START_REF] Bertsekas | Nonlinear programming[END_REF][START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF] that sequentially updates the abundances, the noise variances and the residual term. The proposed Bayesian model and estimation algorithm are validated using synthetic and real hyperspectral images. The obtained results are very promising and show the potential of the proposed mixture and Bayesian models and their associated inference algorithm.

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

MIXTURE MODEL

The LMM is widely used because of its simplicity. However, there are a lot of situations where the linear model is not valid because of the presence of nonlinearity or other mismodelling effects. These effects require the elaboration of a more sophisticated model to improve the abundance estimates. This paper proposes a new formulation based on a residual component model [START_REF] Kalaitzis | Residual components analysis[END_REF] that is expressed as the sum of a linear model and a residual term. The general observation model for the (L × 1) pixel spectrum y i,j is given by

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

Σ = diag σ 2 , σ 2 = σ 2 1 , • • • , σ 2 L
T is an (L × 1) vector containing the noise variances, L is the number of spectral bands, M is a fixed endmember matrix that is assumed known (extracted using an endmember extraction algorithm) and φ N L i,j is a residual term that account for NL interactions. Similarly to [START_REF] Veganzones | A New Extended Linear Mixing Model to Address Spectral Variability[END_REF][START_REF] Shaw | Spectral imaging for remote sensing[END_REF], model [START_REF] Hapke | Bidirectional reflectance spectroscopy. I. Theory[END_REF] includes the effect of illumination variation via the pixel dependent illumination coefficient ci,j (we consider the same parameter for all the endmembers). Moreover, the effect of bilinear interactions between the endmembers is introduced by considering the following residual term

φ N L i,j (M ) = c 2 i,j Q (M ) γ i,j (2)
where

γ i,j = γ (1) i,j , • • • , γ (R) i,j , γ (1,2) i,j , • • • , γ (R-1,R) i,j T , ∀i, j is the (D × 1) vector of positive nonlinear coefficients, D = R(R+1) 2 , Q (M ) = m1 m1, • • • , mR mR, √ 2m1 m2, • • • , √ 2 mR-1 mR)
, is the (L × D) matrix gathering the interaction spectra, and denotes the Hadamard (termwise) product. Note that model ( 1) reduces to the model [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF] for ci,j = 1, ∀, i, j and to the LMM for γ i,j = 0, and ci,j = 1, ∀, i, j. In addition, model ( 1) has a polynomial-like form as for the bilinear models (GBM [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], PP-NMM [START_REF] Altmann | Supervised nonlinear spectral unmixing using a postnonlinear mixing model for hyperspectral imagery[END_REF], Nascimento [START_REF] Bioucas-Dias | Nonlinear mixture model for hyperspectral unmixing[END_REF], Fan [START_REF] Fan | Comparative study between a new nonlinear model and common linear model for analysing laboratory simulated-forest hyperspectral data[END_REF] and Meganem [START_REF] Meganem | Linear-quadratic mixing model for reflectances in urban environments[END_REF] models). Note finally that model (1) (with no illumination variation) has been studied in [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF] when considering a Markov chain Monte-Carlo (MCMC) approach and have shown good performance for processing hyperspectral images. However, the MCMC estimation algorithm was computationally expensive, and we consider in this paper a faster algorithm based on a coordinate descent algorithm.

HIERARCHICAL BAYESIAN MODEL

The unknown parameters associated with model ( 1) are Θ = (A, c, Γ, σ 2 , , where A (resp. Γ) gathers all the abundance vectors (resp. nonlinear coefficients), and is a hyperparameter described later in the text. A hierarchical Bayesian model is considered for the estimation of Θ. This approach is well suited for ill-posed inverse problems, since it accounts for the known parameter constraints using a priori distributions. The next sections present the likelihood and the proposed parameter prior distributions.

Likelihood

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

y i,j |ai,j, ci,j, γ i,j , Σ ∼ N µ i,j , Σ . (3)
where µ i,j = ci,jM ai,j + φ N L i,j (M ) and ∼ means "is distributed according to". The independence between the observed pixels leads to f (Y |Θ) = i j f (y i,j |ai,j, ci,j, γ i,j , Σ), where Y gather all the pixels.

Parameter priors

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

Abundance matrix A

The abundances should satisfy the positivity and sum-to-one (PSTO) constraints. Akin to [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF][START_REF] Dobigeon | Semisupervised linear spectral unmixing using a hierarchical Bayesian model for hyperspectral imagery[END_REF], and since there is no additional information about A, we assign the abundances a uniform prior in the simplex S = ai,j ar,i,j ≥ 0, ∀r and R r=1 ar,i,j = 1 .

Prior for c

In absence of illumination variability, the parameter ci,j has a fixed value #1 for all pixels [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF][START_REF] Heinz | Fully constrained least-squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF][START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF]. This paper accounts for the illumination variability by allowing this parameter to fluctuate around #1. This is achieved by considering the following conjugate Gaussian prior for ci,j ∼ N 1, η 2 , where η 2 is a small fixed variance (η 2 = 0.01 in the rest of the paper). For simplicity, we denote "x|θ ∼ ...", by "x ∼ ..." when the parameter θ is a user fixed parameter. Note finally that the joint prior of c is obtained by assuming a priori independence between the coefficients ci,j, as follows f (c) = i,j f (ci,j).

Nonlinear coefficients γ i,j

Due to physical constraints, the nonlinear coefficients should satisfy the positivity constraint. Similarly to [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF], γ i,j are assigned the following truncated Gaussian prior

γ i,j | 2 i,j ∼ N (R+) D 0D, 2 i,j ID , (4)
where ID denotes the D × D identity matrix and 2 i,j is a variance parameter that is pixel dependent. From (4), it is clear that this variance is related to the strength of the nonlinearities at the pixel (i, j) (via the norm ||γ i,j || 2 ). Moreover, as in [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF], we assume the nonlinear energies to vary smoothly from one pixel to another which will be introduced by considering a specific prior for 2 i,j , as explained in Section 3.3. Note finally that the joint prior of Γ is obtained by assuming a priori independence between the nonlinear coefficients, as follows f (Γ| ) = i,j f γ i,j | 2 i,j .

Noise variances

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

Hyperparameter priors

Due to the spatial organization of hyperspectral images, we expect the energies of the nonlinear coefficients γ i,j to vary smoothly from one pixel to another. This behavior is obtained by introducing an auxiliary variable w (of size Nrow × Ncol) and assigning a gamma Markov random field (GMRF) prior for ( , w) (see [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF][START_REF] Dikmen | Gamma markov random fields for audio source modeling[END_REF] for more details regarding this prior). An interesting property of this joint prior is that the conditional prior distributions of and w reduce to conjugate inverse gamma (IG) and gamma (G) distributions as follows

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

ρ1,i,j(w) = (w 2 i,j + w 2 i+1,j + w 2 i,j+1 + w 2 i+1,j+1 )/4, ρ2,i,j( ) = ( -2 i,j + -2 i-1,j + -2 i,j-1 + -2 i-1,j-1 )/4, (6)
and ζ is a fixed coupling parameter that controls the amount of spatial smoothness enforced by the GMRF. Note that (5) ensures that i+1,j

2 i+1,j+1 2
i+1,j-1 2 i,j 2 i-1,j w 2 i+1,j+1 w 2 i+1,j 2 i,j+1 2 i,j-1 w 2 i,j w 2 i,j+1 2 i-1,j+1 2 i-1,j-1 ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ
Fig. 1. Considered 1st order GMRF neighborhood structure. each 2 i,j is connected to four neighbor elements of w and vice-versa (see Fig. 1). Note also that the energies 2 i,j are conditionally independent and the 1st order neighbors (i.e., the spatial correlation) is introduced via the interaction between and the auxiliary variables w.

Posterior distributions

The joint posterior distribution of the hierarchical Bayesian model can be computed from the following Bayes' rule [START_REF] Bioucas-Dias | Nonlinear mixture model for hyperspectral unmixing[END_REF] where ∝ means "proportional to" and we have assumed a priori independence between the parameters of each model. The MMSE and MAP estimators associated with the posterior [START_REF] Bioucas-Dias | Nonlinear mixture model for hyperspectral unmixing[END_REF] are not easy to determine. In this paper, and akin to [START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF], we propose to evaluate the MAP estimator by using an optimization technique maximizing the posterior (7) w.r.t. the parameters of interest.

f (Θ|Y ) ∝ f (Y |Θ)f (A) f σ 2 f (c) f (γ| ) f ( , w) ,

COORDINATE DESCENT ALGORITHM

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

• Conditional of A: truncated Gaussian distribution (whose one maximum is obtained with SUNSAL-FCLS 1 [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF])

• Conditional of Γ: positive truncated Gaussian distribution (whose one maximum is obtained with SUNSAL-CLS [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF])

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

1 SUNSAL-FCLS satisfies the PSTO constraints while SUNSAL-CLS only ensure the positivity constraint.

• Conditional of w : gamma distribution (analytical expression of the maximum) • Conditional of σ 2 : inverse gamma distribution (analytical expression of the maximum) • Conditional of c: Gaussian distribution (analytical expression of the mean).

Regarding the sequence generated by the coordinate descent algorithm, the proposition 2.7.1 in [START_REF] Bertsekas | Nonlinear programming[END_REF] asserts that its limit points are stationary points of [START_REF] Bioucas-Dias | Nonlinear mixture model for hyperspectral unmixing[END_REF] provided that the maximum of that function w.r.t. Θ along each coordinate is unique. This is easily checked for all the parameters except for c NL . Indeed, the cost function writes as a 4-order polynomial w.r.t. c NL (leading to 3 possible maxima) and we have chosen the one that maximizes it in the interval [0.2, 3]. Note also that the cost function is not convex, thus, the solution obtained might depend on the initial values that need to be chosen carefully. Therefore, the abundances A are initialized with SUNSAL-FCLS [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF], the nonlinear coefficients are initialized by 0, the noise variance is initialized by HySime [START_REF] Bioucas-Dias | Hyperspectral subspace identification[END_REF], the illumination coefficient c is initialized by considering the sum of the abundances that are estimated using only the positivity constraint with SUNSAL-CLS [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF]. With these initializations, the proposed algorithm reached minima of good quality in the considered simulations (see Sections 5 and 6).

SIMULATION RESULTS ON SYNTHETIC DATA

This section evaluates the performance of the proposed CDA-NL algorithm on a synthetic image denoted by I. This image contains 100 × 100 pixels, L = 207 spectral bands and has been generated using R = 3 endmembers extracted from the ENVI software library. The image has also been corrupted by i.i.d. Gaussian noise for a fair comparison with SU algorithms using this assumption. To investigate the behavior of the proposed approach w.r.t. different mixture models, image I has been generated according to 4 linear/nonlinear models. For this, an image partition into 4 classes has been generated by considering a Potts-Markov random field. The four spatial classes are associated with the LMM, NL model (1) (with 2 = 0.1), GBM (with random nonlinear coefficients in [0.8; 1]) and PPNMM (with b = 0.5), respectively. Note that the abundances have been clustered inside the simplex using a Dirichlet distribution and that the illumination coefficient c increases linearly from the left to the right of the image in [0.9, 1.15]. Table 1 compares the proposed CDA-NL algorithm with the FCLS [START_REF] Heinz | Fully constrained least-squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF], SUNSAL-CLS [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF], SKhype [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinearfluctuation model[END_REF], and RCA-MCMC [START_REF] Altmann | Residual component analysis of hyperspectral images: Application to joint nonlinear unmixing and nonlinearity detection[END_REF] algorithms. The criteria used are the quadratic distance between the actual and estimated abundances (RMSE), the reconstruction error (RE) and the spectral angle mapper (SAM) between the observed and reconstructed pixel spectra. Table 1 highlights the good performance of CDA-NL for the different mixture models and shows that CDA-NL has the best overall results. It also shows that CDA-NL presents a good computational times that is competitive with the NL unmixing algorithms SKhype and RCA-MCMC.

RESULTS ON REAL DATA

This section illustrates the performance of the proposed algorithm when applied to a real hyperspectral image. This image was acquired over Moffett Field, CA, in 1997 by the AVIRIS. The considered dataset contains 100 × 100 pixels, L = 152 spectral bands (after removing water absorption bands) acquired in the interval 0.4 -2.5µm, has a spatial resolution of 100m and is mainly composed of three components: water, soil, and vegetation (see Fig. 2 (Top-left)). The VCA algorithm [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF] was used to estimate R = 3 endmembers. The estimated CDA-NL abundance maps are in good agreement with the state-of-the-art algorithms and are not shown here for brevity (see [START_REF] Halimi | Hyperspectral unmixing in presence of endmember variability, nonlinearity or mismodelling effects[END_REF] for more details). .

CONCLUSIONS

This paper introduced a hyperspectral mixture model and its associated Bayesian algorithm for nonlinear hyperspectral unmixing. The proposed model considered the spatial illumination variation and the presence of nonlinearity effects. A hierarchical Bayesian model was proposed to introduce the known constraints on the parameters of interest. Those parameters were estimated using a coordinate descent algorithm that showed a reduced computational cost when compared to state-of-the-art algorithms. The proposed algorithm showed good performance when processing synthetic data generated with the linear and nonlinear models. Results on real data confirmed the good performance of the proposed algorithm and showed its ability to extract different features in the observed scene. Future work includes the estimation of the hyperparameters associated with the proposed Bayesian model.

Fig. 2 (

 2 top-right) shows the energies of the difference between the reconstructed signal and the linear model. It is clear that the residual components are mainly located in the coastal region and in water. Figs. 2 (bottom-left) and (bottom-right) show the effects of the illumination variation and the nonlinear term, respectively. It can be seen that the illumination variation mainly occurs in water while the nonlinear effects are located in the coastal zone and in presence of vegetation.

Fig. 2 .

 2 Fig. 2. (Top-left) Real Moffett image. (Top-right) Square root of the energies of the difference between the reconstructed signal and the linear model (||ŷ i,j -M âi,j||). (Bottom-left) estimated illumination variation |1 -ci,j|. (Bottom-right) square root of the energies of the residual terms φ N L i,j

Table 1 .

 1 Results on synthetic data.

		RMSE (classes) (×10 -2 ) C2 C3 LMM RCA-NL GBM PPNMM C1 C4	RMSE (×10 -2 ) (×10 -2 ) (×10 -2 ) RE SAM	Time (s)
	FCLS	10.24	44.72	15.48	23.98	24.76	15.74	10.64	1.7
	SUNSAL-CLS	3.84	33.81	5.68	8.45	16.55	4.17	7.57	0.07
	SKhype	1.67	11.92	2.21	2.81	5.87	-	-	547
	RCA-MCMC	5.87	6.29	5.44	3.93	5.66	-	-	9009
	CDA-NL	1.62	7.27	2.16	2.89	3.86	2.86	6.16	430

This work was supported in part by the EPSRC Grants EP/JO15180/1, EP/N003446/1, and EP/K015338/1, and in part by the HYPANEMA ANR Project under Grant ANR-12-BS03-003.