Improved Hyperspectral Unmixing With Endmember Variability Parametrized Using an Interpolated Scaling Tensor

TL;DR

This paper introduces a two-step hyperspectral unmixing method that leverages pure pixel information and low-rank spectral variability modeling to improve unmixing accuracy and robustness.

Contribution

It proposes a novel approach that separates the estimation of spectral variability from abundance calculation, enhancing stability and exploiting pure pixel information effectively.

Findings

Improved unmixing accuracy on synthetic and real data.

Effective modeling of spectral variability as a smooth, low-rank tensor.

Enhanced robustness against ill-posed optimization problems.

Abstract

Endmember (EM) variability has an important impact on the performance of hyperspectral image (HI) analysis algorithms. Recently, extended linear mixing models have been proposed to account for EM variability in the spectral unmixing (SU) problem. The direct use of these models has led to severely ill-posed optimization problems. Different regularization strategies have been considered to deal with this issue, but none so far has consistently exploited the information provided by the existence of multiple pure pixels often present in HIs. In this work, we propose to break the SU problem into a sequence of two problems. First, we use pure pixel information to estimate an interpolated tensor of scaling factors representing spectral variability. This is done by considering the spectral variability to be a smooth function over the HI and confining the energy of the scaling tensor to a low-rank structure. Afterwards, we solve a matrix-factorization problem to estimate the fractional abundances using the variability scaling factors estimated in the previous step, what leads to a significantly more well-posed problem. Simulation swith synthetic and real data attest the effectiveness of the proposed strategy.

Tables

Table 1. Table 1 : Simulations with synthetic and real data.

Synthetic HI					Houston HI
	${\text{RMSE}}_{𝑨}$	${\text{RMSE}}_M$	${\text{SAM}}_M$	${\text{RMSE}}_{𝑹}$	${\text{RMSE}}_{𝑹}$
FCLS	0.0416	–	–	0.0212	0.0478
PLMM	0.0353	0.0220	0.0232	0.0117	0.0360
ELMM	0.0373	0.0135	0.0174	0.0121	0.0033
GLMM	0.0343	0.0128	0.0163	0.0118	0.0003
Proposed	0.0233	0.0130	0.0164	0.0123	0.0085

Equations

$${\mathbf{r}}_n=\mathbf{M}{\mathbf{\alpha }}_n+{\mathbf{e}}_n,\text{subject to }{1}^{\top }{\mathbf{\alpha }}_n=1\text{ and }{\mathbf{\alpha }}_n⪰0$$

$${\mathbf{r}}_n=(\mathbf{M}\odot {\mathbf{\Psi }}_n){\mathbf{\alpha }}_n+{\mathbf{e}}_n$$

$$J(\mathbf{A},\text{amsmathbb}\mathbf{M})=\frac{1}{2}\Vert \mathbf{R}-f(\mathbf{A},\text{amsmathbb}\mathbf{M}){\Vert }_F^2+\mathcal{R}(\mathbf{A})+\mathcal{R}(\text{amsmathbb}M)$$

$$J_{\Psi }(\Psi ,\Phi )$$

$$+{\lambda }_{\Psi }\sum _{k=1}^R\sum _{(n_1,n_2)\in {\mathcal{N}}_{\mathcal{P},k}}\Vert {\tilde{\mathbf{r}}}_{n_1,n_2}-{\mathbf{m}}_{0,k}\odot [\Psi {]}_{n_1,n_2,:,k}{\Vert }^2$$

$$\text{subject to }\mathrm{rank}(\Phi )=r$$

$$\underset{\Psi ,\Phi }{\min }{J}_{\Psi }(\Psi ,\Phi )$$

$$\underset{\Psi }{\min }$$

$$+{\lambda }_{\Psi }\sum _{k=1}^R\sum _{(n_1,n_2)\in {\mathcal{N}}_{\mathcal{P},k}}\Vert {\tilde{\mathbf{r}}}_{n_1,n_2}-{\mathbf{m}}_{0,k}\odot [\Psi {]}_{n_1,n_2,:,k}{\Vert }^2$$

$$[\hat{\Psi }{]}_{n_1,n_2,\ell ,k}$$

$$=\{\begin{array}{ll}\frac{[\Phi {]}_{n_1,n_2,\ell ,k}+\epsilon }{1+\epsilon },& (n_1,n_2)\notin {\mathcal{N}}_{\mathcal{P},k}\\ \frac{[\Phi {]}_{n_1,n_2,\ell ,k}+\epsilon +{\lambda }_{\Psi }{m}_{\ell ,k}^0{\tilde{r}}_{\ell ,n_1,n_2}}{1+\epsilon +{\lambda }_{\Psi }(m_{\ell ,k}^0{)}^2},& (n_1,n_2)\in {\mathcal{N}}_{\mathcal{P},k}\end{array}$$

$$\underset{\Phi }{\min }\Vert \Psi -\Phi {\Vert }_F^2\text{subject to }\mathrm{rank}(\Phi )=r$$

$$\Phi =\sum _{i=1}^r{\xi }_i{\mathbf{z}}_i^{(1)}\circ {\mathbf{z}}_i^{(2)}\circ {\mathbf{z}}_i^{(3)}\circ {\mathbf{z}}_i^{(4)}.$$

$$\displaystyle\Big{(}\widehat{{\boldsymbol{\Xi}}}\,\,,$$

$$\displaystyle\mathop{\arg\min}_{{\boldsymbol{\Xi}},\boldsymbol{Z}^{(1)},\boldsymbol{Z}^{(2)},\boldsymbol{Z}^{(3)},\boldsymbol{Z}^{(4)}}\Big{\|}\mathbb{\Psi}-\sum_{i=1}^{r}\xi_{i}\boldsymbol{z}_{i}^{(1)}\circ\boldsymbol{z}_{i}^{(2)}\circ\boldsymbol{z}_{i}^{(3)}\circ\boldsymbol{z}_{i}^{(4)}\Big{\|}^{2}_{F}$$

$$\widehat{\mathbb{\Phi}}=\big{\ldbrack}\widehat{{\boldsymbol{\Xi}}}\,\,;\,\widehat{\boldsymbol{Z}}^{(1)};\widehat{\boldsymbol{Z}}^{(2)};\widehat{\boldsymbol{Z}}^{(3)}\big{\rdbrack}.$$

$$J(\mathbf{A},\text{amsmathbb}\mathbf{M})=$$

$$\displaystyle+\lambda_{A}\big{(}\|{\mathcal{H}}_{h}(\boldsymbol{A})\|_{2,1}+\|{\mathcal{H}}_{v}(\boldsymbol{A})\|_{2,1}\big{)}$$

$$\text{subject to }\mathbf{A}⪰0,{\mathbf{A}}^{\top }1=1,\text{amsmathbb}\mathbf{M}⪰0$$

$$\displaystyle\mathop{\min}_{\amsmathbb{M}\succeq{\boldsymbol{0}}}\,\,\frac{1}{2}\!\sum_{n=1}^{N}\!\big{(}\|\boldsymbol{r}_{n}-\boldsymbol{M}_{n}\boldsymbol{\alpha}_{n}\|^{2}\!+\!\lambda_{M}\|\boldsymbol{M}_{n}-\boldsymbol{M}_{0}\odot\boldsymbol{\Psi}_{n}\|^{2}_{F}\big{)}.$$

$${\hat{\mathbf{M}}}_n=({\mathbf{r}}_n{\mathbf{\alpha }}^{\top }+{\lambda }_M{\mathbf{M}}_0\odot {\mathbf{\Psi }}_n)({\mathbf{\alpha }}_n{\mathbf{\alpha }}_n^{\top }+{\lambda }_M{\mathbf{I}}_R{)}^{-1}$$

$$\displaystyle\mathop{\min}_{\boldsymbol{A}}\,\,\frac{1}{2}\!\sum_{n=1}^{N}\|\boldsymbol{r}_{n}-\boldsymbol{M}_{n}\boldsymbol{\alpha}_{n}\|^{2}\!+\lambda_{A}\big{(}\|{\mathcal{H}}_{h}(\boldsymbol{A})\|_{2,1}+\|{\mathcal{H}}_{v}(\boldsymbol{A})\|_{2,1}\big{)}$$

$$\text{subject to }\mathbf{A}⪰0,{\mathbf{A}}^{\top }1=1$$

$${\text{RMSE}}_{\text{amsmathbb}X}=\sqrt{\frac{1}{N_{\text{amsmathbb}X}}\Vert \text{amsmathbb}X-{\text{amsmathbb}X}^{\ast }{\Vert }_F^2}$$

$${\text{SAM}_{\amsmathbb{M}}=\frac{1}{N}\sum_{n=1}^{N}\sum_{k=1}^{R}\arccos\bigg{(}\frac{\boldsymbol{m}_{k,n}^{\top}\boldsymbol{m}_{k,n}^{*}}{\|\boldsymbol{m}_{k,n}\|\|\boldsymbol{m}_{k,n}^{*}\|}\bigg{)}.}$$

Full text

**Index Terms— ** Hyperspectral data, endmember variability, GLMM, pure pixel, tensor interpolation.

1 Introduction

Spectral Unmixing (SU) is part of a number of algorithms that retrieve vital information from hyperspectral images (HIs) in many applications [1]. SU aims at extracting the spectral signatures of materials present in the HI of a scene, as well as the proportion in which they contribute to each HI pixel. Many parametric models have been proposed to describe the interaction between light and the target surface [1, 2]. The simplest of such models is the Linear Mixing Model (LMM), which considers that the observed reflectance of an HI pixel is obtained from a convex combination of the spectral signatures of pure materials. This model imposes a convex geometry to the SU problem, where all HI pixels are confined to a simplex whose vertices are the reflectances of pure materials, usually termed endmembers (EMs). The linearity and convexity of the LMM model lead to an interpretation of its coefficients as the relative abundances of each pure material in the HI. Nevertheless, some characteristics of practical HIs cannot be modeled by the standard LMM, such as nonlinearities [2, 3, 4, 5] or variations of the EMs along the image [6, 7, 8]. More sophisticated models are required when such nonidealities have important impact on the HI.

EM variability can result, for instance, from environmental conditions, illumination, atmospheric or temporal changes [9]. Its occurrence may incur the propagation of significant estimation errors throughout the unmixing process [6]. Most of the methods proposed so far to deal with spectral variability can be classified in three major groups: endmembers as sets, endmembers as statistical distributions and, more recently, methods that incorporate the variability in the mixing model, often using physically motivated concepts [10]. The method proposed in this work combines elements from the first and third groups. Specifically, we adopt a modified mixing model over a very specific variational set.

In the first group, EM variability is addressed by considering different types of variational sets, among which spectral bundles are prominent [9]. One nice property of these methods is that the bundles can be directly extracted from the observed HI. Few works expand this idea by adopting pixel-dependent EMs obtained using some kind of spatial interpolation of the original EM signatures [11, 12, 13]. These interpolated EMs sets are then used to unmix the data using the standard LMM. This type of approach often profits from a priori information about EM spectra (pure pixels) and their position in the image. Interpolation strategies include linear regression [11, 14, 15] and kriging [12, 13]. However, these strategies lack flexibility to adapt to variations in the endmember signatures that lie beyond the descriptive capability of the spectral bundles.

In the third group, different extensions of the LMM have been proposed to cope with spectral variability [6, 7, 8]. The models differ with respect to the type of variability they represent and their physical motivation. The Perturbed LMM model (PLMM) [6] introduces an additive perturbation to the EM matrix. It allows the modeling of arbitrary spectral variations, but lacks physical motivation. The Extended LMM (ELMM) [7] and its generalization, the Generalized LMM (GLMM) [8], are physically motivated, but differ in their ability to model arbitrary variability. The ELMM performs well if spectral variability is due mainly to illumination variations, but it lacks flexibility when the EMs are subject to more complex spectral distortions. The GLMM is able to model arbitrary EM variability by considering band-dependent multiplicative factors, thus linking the amount of spectral variability to the magnitude of the EM reflectance in each band. This effect is consistent with empirical observations, leading to a more physically reasonable model. One drawback of extended LMMs is that they lead to ill-posed optimization problems. Hence, regularization strategies must be devised to yield meaningful results [16, 17]. Different regularization strategies have been recently proposed. For instance, [18] and [19] proposed a data-dependent multiscale strategy to exploit the spatial uniformity existing in HIs. The approaches in [20, 21] assume that most of the energy of HIs is confined to a low-dimensional structure within the image space, and exploit low-rank tensor decomposition techniques to impose regularity to the solution.

In this work we propose a hybrid approach by leveraging pure pixel information in the context of parametric endmember models. Specifically, we propose to break the SU problem into two sequential problems. We first use pure pixel information to estimate an interpolated tensor of spectral variability scaling factors by assuming that its energy is confined to a low-rank structure. Afterwards, we solve a matrix-factorization problem to estimate the fractional abundances. This is done by using the variability scaling factors estimated in the previous step in the form of a regularization term, which leads to a significantly more well-posed problem. We adopt the GLMM to parametrically model spectral variability, as it combines flexibility and physical motivation.

This paper is organized as follows. Section 2 briefly reviews the LMM and its GLMM extended version. Section 3 introduces the proposed formulation and the corresponding SU algorithm. The performance of the proposed method is compared with competing algorithms in Section 4. Finally, the conclusions are presented in Section 5.

2 Linear Mixing Models

The Linear Mixing Model (LMM) [1] assumes that a given a pixel $\boldsymbol{r}_{n}=[r_{n,1},\,\ldots,\,r_{n,L}]^{\top}$, with $L$ bands, is represented as

[TABLE]

where $\boldsymbol{M}\in\amsmathbb{R}^{L\times R}$ is a matrix whose columns are the $R$ endmember spectral signatures $\boldsymbol{m}_{k}$, $\boldsymbol{\alpha}_{n}$ is the abundance vector, $\boldsymbol{e}_{n}$ is an additive white Gaussian noise (WGN) and $\succeq$ is the entrywise $\geq$ operator. The LMM assumes that the endmember spectra are fixed for all pixels $\boldsymbol{r}_{n}$, $n=1,\ldots,N$, in the HI. This assumption jeopardizes the accuracy of estimated abundances in many circumstances due to the spectral variability existing in a typical scene.

To mitigate this issue, the GLMM [8] models arbitrary EM variability by considering band-dependent multiplicative factors, which links the amount of spectral variability to the magnitude of the EM reflectance in each band.

The GLMM introduces a scaling matrix ${\boldsymbol{\Psi}}_{n}\in{\amsmathbb{R}^{L\times R}}$ with entries $[{\boldsymbol{\Psi}}_{n}]_{\ell,k}=\psi_{n_{\ell,k}}\geq 0$, which acts individually on each wavelength. The model represents the $n$-th HI pixel as

[TABLE]

where $\odot$ stands for the Hadamard product.

3 Proposed Unmixing Algorithm

Given a parametric extended linear mixing model $f$, and an HI $\boldsymbol{R}=[\boldsymbol{r}_{1},\ldots,\boldsymbol{r}_{N}]$, the SU problem can be generally cast as the determination of $\boldsymbol{A}$ and $\amsmathbb{M}$ that minimize a risk functional of the form

[TABLE]

where $\amsmathbb{M}$ is a 3-D tensor obtained by stacking all pixel-dependent endmember matrices $\boldsymbol{M}_{n}$, such that $[\amsmathbb{M}]_{n,:,:}=\boldsymbol{M}_{n}$, $\boldsymbol{A}$ is the abundance matrix, and $\mathcal{R}(\boldsymbol{A})$ and $\mathcal{R}(\amsmathbb{M})$ are regularization terms to improve the problem conditioning. Although parametric EM models such as the ELMM and GLMM have enough flexibility to adequately represent the variability in the endmember spectra, the resulting unmixing problem is highly ill-posed. This problem is accentuated by the fact that previous formulations such as [16, 8] only considered a single reference, or average endmember matrix, as a priori information (or initialization), while estimating the variability maps and endmembers for all pixels in the scene. Although good performance has been reported in [16, 8], the correct estimation of the parametric models for complex scenes remains a challenge.

A characteristic of many scenes, which remains unexplored in this context, is the existence of multiple pure pixels in an observed HI. This information can be leveraged to help in the estimation of the parametric endmember model, introducing more information into the problem and reducing its ill-posedness.

We propose to break the unmixing problem into a sequence of two problems:

• (i)

Using pure pixel information extracted from the HI by standard EM extraction algorithms and assuming smooth spectral variability throughout the HI, estimate an interpolated tensor of variability factors $\mathbb{\Psi}$ whose energy is assumed to be confined to a low-rank structure.

• (ii)

Using the estimated tensor as a prior information in a standard matrix-factorization problem, solve it to estimate the endmember and abundance matrices.

3.1 Estimating the Tensor of Variability Factors

The objective of this first problem is to estimate the scaling matrices $\boldsymbol{\Psi}_{n}$, $n=1,\ldots,N$, in (2). To this end, we assume:

• a)

the availability of a reference endmember matrix $\boldsymbol{M}_{0}$, which can be obtained using any endmember extraction method;

• b)

the knowledge of the set $\mathcal{N}_{\mathcal{P},k}$ of locations $(n_{1},n_{2})$ of pure pixels for the $k$-th endmember, for all $k=1,\dots,R$; 111Pure pixels are here defined here as a set of pixels whose spectral distance relative to the reference EMs in $\boldsymbol{M}_{0}$ is less than a specified threshold.

• c)

that the EM variability changes smoothly in the HI, so that each endmember variability function has most of its energy confined to a low-dimensional structure in the image space.

Consider the 4-D tensor $\mathbb{\Psi}\in\amsmathbb{R}^{N_{1}\times N2\times L\times R}$, $N_{1}\times N_{2}=N$, of variability weights with entries $[\mathbb{\Psi}]_{n_{1},n_{2},\ell,k}=\psi_{n_{1},n_{2},\ell,k}$, corresponding to the scaling factor to be applied to the $\ell$-th band of the $k$-th endmember in the pixel at position $(n_{1},n_{2})$ of the HI. $\mathbb{\Psi}$ generalizes the GLMM scaling factors. We propose to estimate the variability tensor $\mathbb{\Psi}$ by solving the following optimization problem:

[TABLE]

In (4), $\mathbb{\Phi}$ is a low-rank tensor which imposes a low-rank constraint on the variability maps to enforce the smooth variation. $\mathbb{1}$ is an $N_{1}\times N_{2}\times L\times R$ tensor of ones, so that $\|\mathbb{\Psi}-\mathbb{1}\|^{2}_{F}$ is a regularization term that enforces the variability to be small, and whose importance is controlled by the parameter $\varepsilon>0$. $\widetilde{\boldsymbol{r}}_{n_{1},n_{2}}$ represents the pixel at position $(n_{1},n_{2})$, $\boldsymbol{m}_{0,k}$ is the $k$-th column of $\boldsymbol{M}_{0}$ with elements $m^{0}_{\ell,k}$, and the third term enforces the similarity between pure pixels at known positions $(n_{1},n_{2})$ and their corresponding endmembers and variability maps. Its strictness is controlled by $\lambda_{\mathbb{\Psi}}>0$. The associated optimization problem is given by

[TABLE]

This problem is non-convex and generally NP-hard due to the rank constraint. We propose to solve it using the alternating optimization approach described in Algorithm 1, as follows.

3.1.1 Optimizing with respect to $\mathbb{\Psi}$

To solve problem (5) with respect to $\mathbb{\Psi}$ we consider $\mathbb{\Phi}$ fixed. Thus, the optimization problem with respect to $\mathbb{\Psi}$ is given by

[TABLE]

By taking the derivative with respect to each position $[\mathbb{\Psi}]_{n_{1},n_{2},\ell,k}$ of $\mathbb{\Psi}$ and setting it equal to zero, we end up with the following solution:

[TABLE]

for $\ell=1,\ldots,L$, $k=1,\ldots,R$ and $n_{i}=1,\ldots,N_{i}$,   $i=1,2$.

3.1.2 Optimizing with respect to $\mathbb{\Phi}$

The optimization problem with respect to $\mathbb{\Phi}$ is given by

[TABLE]

where we write tensor $\mathbb{\Phi}$ as

[TABLE]

where $\circ$ denotes the outer product and $r$ is a small number since we assume that most of the energy of $\mathbb{\Psi}$ lies in a low-rank structure. Using (9) in (8) leads to the following optimization problem:

[TABLE]

where {\boldsymbol{\Xi}}=\text{Diag}_{4}\big{(}\xi_{1},\ldots,\xi_{r}\big{)} is an order-4 diagonal tensor with $[{\boldsymbol{\Xi}}]_{i,i,i,i}=\xi_{i}$. Problem (10) can be solved using an alternating least-squares strategy [22]. The solution $\widehat{\mathbb{\Phi}}$ is then obtained using the full multilinear product [21]:

[TABLE]

3.2 The Spectral Unmixing Problem

Given a reference EM matrix $\boldsymbol{M}_{0}$ and the EM scaling factors $\boldsymbol{\Psi}_{n}=\widehat{\mathbb{\Psi}}_{n,:,:}$ for each pixel (with each $n$ corresponding to a single pair $(n_{1},n_{2})$), the SU problem can be formulated as the minimization of the risk functional (3) with respect to $\boldsymbol{A}$ and $\amsmathbb{M}$, using an appropriate EM model and regularization terms. We propose to minimize the following regularized cost functional:

[TABLE]

The last term in (3.2) promotes spatial regularity in the abundance maps. $\mathcal{H}_{h}$ and $\mathcal{H}_{v}$ are linear operators that compute the horizontal and vertical gradients of a bidimensional signal, acting separately for each material of $\boldsymbol{A}$, and $\|\boldsymbol{X}\|_{2,1}=\sum_{n=1}^{N}\|\boldsymbol{x}_{n}\|_{2}$ is the $\mathcal{L}_{2,1}$ norm.

Cost function in (3.2) is non-smooth and non-convex with respect to both $\boldsymbol{A}$ and $\amsmathbb{M}$, but convex with respect to each of them. Following the approaches used in [7, 8], we search for a locally optimal solution by alternately minimizing (3.2) with respect to each variable, as described in Algorithm 2.

3.2.1 Optimization with respect to $\amsmathbb{M}$

Considering only the terms in (3.2) that depend on $\amsmathbb{M}$, the optimization problem for the EM tensor becomes

[TABLE]

Problem (13) can be solved individually for each pixel $\boldsymbol{r}_{n}$. Relaxing the nonnegativity constraint on $\amsmathbb{M}$, the solution is given by

[TABLE]

An approximate solution to the original constrained problem (13) can then be obtained by projecting $\,\widehat{\!\boldsymbol{M}}_{n}$ onto the nonnegative orthant $\amsmathbb{R}_{+}^{L\times R}$ by attributing zero to any negative entries [7, 8].

3.2.2 Optimization with respect to $\boldsymbol{A}$

Considering only the terms in (3.2) that depend on $\boldsymbol{A}$, the optimization problem for the abundance matrix becomes

[TABLE]

This problem is non-smooth and not separable with respect to the different pixels in the HI. Nevertheless, problem (3.2.2) can be solved efficiently using the Alternating Direction Method of the Multipliers (ADMM) [23]. The procedure is described in detail in [7].

4 Experimental Results

We now illustrate the performance of the proposed method through simulations with both synthetic and real data. We compare the proposed method with the fully constrained least squares (FCLS), the ELMM [7], the PLMM [6] and the GLMM [8]. In all experiments, the reference EM matrix $\boldsymbol{M}_{0}$ was extracted from the observed HI using the VCA algorithm [24]. The sets $\mathcal{N}_{\mathcal{P},k}$ of pure pixel locations were estimated by selecting the pixels with the smallest spectral angle relative to the reference EMs in $\boldsymbol{M}_{0}$. The performances were evaluated using the Root Means Squared Error (RMSE) between the estimated abundance maps ($\text{RMSE}_{\boldsymbol{A}}$), between the EM matrices ($\text{RMSE}_{\amsmathbb{M}}$) and between the reconstructed images ($\text{RMSE}_{\boldsymbol{R}}$). The RMSE between two generic tensors is defined as

[TABLE]

where $N_{\amsmathbb{X}}$ denotes the number of elements in the tensor $\amsmathbb{X}$.

We consider also the Spectral Angle Mapper (SAM) to evaluate the estimated EM tensor

[TABLE]

4.1 Synthetic data

A synthetic dataset was created from spatially correlated abundance maps containing a significant amount of pure pixels. The HI was constructed using reference EMs extracted from the USGS Spectral Library [25]. Spectral variability was introduced following the GLMM model in Eq. (2) to generate pixel-dependent EM signatures with spatial and spectral correlation imposed on $\boldsymbol{\Psi}_{n}$ using a 3-D Gaussian filter. Then, WGN was added to yield a 30dB SNR.

We selected the optimal parameters for each algorithm by performing a grid search using the ranges of parameters suggested by the authors in the original publications. For PLMM we used $\gamma=1$ and searched for $\alpha$ and $\beta$ in the range $\{0.1,\,25\}$ and $\{10^{-9},\,10^{-3}\}$, respectively. For ELMM, GLMM and the proposed method, we selected the parameters in the following ranges: $\lambda_{S},\,\lambda_{M}\in\{0.01,50\}$, $\lambda_{A}\in\{0.001,\,0.1\}$, and $\lambda_{\psi},\,\lambda_{\mathbb{\Psi}}\in\{10^{-4},\,10^{3}\}$. We also used $\varepsilon=10^{-5}$ and $r=10$ fixed in Algorithm 1. Finally, the number of pure pixels extracted from the image was $500$, $100$ and $10$ for the first, second and third EMs, respectively.

The results are shown in Table 1. It can be seen that the proposed method yields a significantly smaller $\text{RMSE}_{\boldsymbol{A}}$ than the other algorithms. It also yields the second-best results $\text{RMSE}_{\amsmathbb{M}}$ and $\text{SAM}_{\amsmathbb{M}}$, both very close to the best results obtained by the GLMM. The reconstruction errors are smaller for PLMM, GLMM and ELMM, which is expected since these methods have a larger number of degrees of freedom, while the proposed method estimates the variability scaling factors a priori. However, smaller reconstruction errors do not necessarily imply better abundance estimation, as evidenced by the $\text{RMSE}_{\boldsymbol{A}}$ results. The execution time of the proposed method is $1.95\times\text{Time}_{\text{GLMM}}$, a modest increase in complexity when compared to the benefits in terms of estimation accuracy.

4.2 Real data

For simulations with real data we considered the Houston dataset, which is known to have four EMs [7, 8]. The parameters for the proposed method were set as $\lambda_{M}=0.1$, $\lambda_{A}=0.01$, $\lambda_{\mathbb{\Psi}}=10^{3}$, $\varepsilon=10^{-5}$ and $r=10$. For the other algorithms we used the same parameters [8]. The sets $\mathcal{N}_{\mathcal{P},k}$ were constructed by extracting all pixels in the HI with an angle smaller than $\cos^{-1}(0.01)$ to the reference EMs. The reconstructed abundance maps for all tested algorithms are shown in Fig. 1, while the reconstruction errors are presented in Table 1.

It can be seen that the proposed method provides abundance estimation that is accurate and comparable with the results obtained using ELMM or GLMM. Furthermore, the proposed method results in stronger concrete and vegetation components in the stadium stands and in the field, respectively, when compared to the other algorithms. As in the synthetic data example, the results in Table 1 show smaller $\text{RMSE}_{\boldsymbol{R}}$ values for GLMM and ELMM. Again, this is not directly related to better abundance estimation. The proposed method demanded a computational time equal to $2.77\times\text{Time}_{\text{GLMM}}$.

5 Conclusions

We have proposed a new spectral unmixing algorithm accounting for spectral variability. Modeling the spectral variability as a smooth function over the image, pure pixel information was leveraged from the HI to estimate a tensor of EM scaling factors prior to unmixing. The use of this predetermined tensor has led to a much more well-posed SU problem when compared to current algorithms, in which the scaling factors are blindly estimated from the image. Simulations using synthetic and real data indicate that the proposed method can lead to significant improvements in abundance estimation accuracy.