Deep Generative Endmember Modeling: An Application to Unsupervised Spectral Unmixing

TL;DR

This paper introduces a deep generative model for spectral endmember variability in hyperspectral images, enabling unsupervised spectral unmixing that adapts to data-driven variability representations.

Contribution

It proposes a novel deep generative endmember model learned via a variational autoencoder, improving flexibility and low-dimensional representation for spectral unmixing.

Findings

Outperforms state-of-the-art algorithms on synthetic data

Effective in modeling spectral variability with a deep generative approach

Unsupervised learning from pure pixel data enhances unmixing accuracy

Abstract

Endmember (EM) spectral variability can greatly impact the performance of standard hyperspectral image analysis algorithms. Extended parametric models have been successfully applied to account for the EM spectral variability. However, these models still lack the compromise between flexibility and low-dimensional representation that is necessary to properly explore the fact that spectral variability is often confined to a low-dimensional manifold in real scenes. In this paper we propose to learn a spectral variability model directly form the observed data, instead of imposing it \emph{a priori}. This is achieved through a deep generative EM model, which is estimated using a variational autoencoder (VAE). The encoder and decoder that compose the generative model are trained using pure pixel information extracted directly from the observed image, what allows for an unsupervised formulation. The proposed EM model is applied to the solution of a spectral unmixing problem, which we cast as an alternating nonlinear least-squares problem that is solved iteratively with respect to the abundances and to the low-dimensional representations of the EMs in the latent space of the deep generative model. Simulations using both synthetic and real data indicate that the proposed strategy can outperform the competing state-of-the-art algorithms.

Tables

Table 1. Table I: Encoder network architecture.

Layer	Activation Function	Number of units
Input	—	$L$
Hidden # 1	ReLU	$\lceil 1.2\times L\rceil +5$
Hidden # 2	ReLU	$\max \{\lceil L/4\rceil ,K+2\}+3$
Hidden # 3	ReLU	$\max \{\lceil L/10\rceil ,K+1\}$

Table 2. Table II: Decoder network architecture.

Layer	Activation Function	Number of units
Hidden # 1	ReLU	$\max \{\lceil L/10\rceil ,K+1\}$
Hidden # 2	ReLU	$\max \{\lceil L/4\rceil ,K+2\}+3$
Hidden # 3	ReLU	$\lceil 1.2\times L\rceil +5$
Output	Sigmoid	$L$

Table 3. Table III: Simulation results using synthetic data.

Data Cube 1 – DC1
	${\text{NRMSE}}_{𝑨}$	${\text{NRMSE}}_{𝕄}$	${\text{SAM}}_{𝕄}$	${\text{NRMSE}}_{𝒀}$	Time [s]
FCLS	0.2854	—	—	0.0350	0.71
PLMM	0.2604	0.1075	0.0440	0.0007	122.09
ELMM	0.2554	0.1032	0.0398	0.0321	8.82
GLMM	0.2480	0.1036	0.0355	0.0235	23.74
DeepGUn	0.0566	0.0944	0.0233	0.0448	75.20
Data Cube 2 – DC2
FCLS	0.1294	—	—	0.0393	0.38
PLMM	0.1197	0.0481	0.0378	0.0336	41.31
ELMM	0.1110	0.0566	0.0382	0.0231	20.25
GLMM	0.1146	0.0534	0.0367	0.0226	17.03
DeepGUn	0.0969	0.0463	0.0323	0.0384	36.40
Data Cube 3 – DC3
FCLS	0.2606	—	—	0.0542	0.34
PLMM	0.2028	0.0928	0.0385	0.0302	59.88
ELMM	0.1997	0.0640	0.0188	0.0238	17.99
GLMM	0.1841	0.0638	0.0185	0.0226	25.74
DeepGUn	0.1613	0.0600	0.0172	0.0457	48.96
Data Cube 4 – DC4
FCLS	0.5109	—	—	0.1712	0.50
PLMM	0.5066	0.6245	0.4874	0.0320	269.48
ELMM	0.4385	0.4712	0.1451	0.0106	18.36
GLMM	0.4371	0.4855	0.1972	0.0108	21.15
DeepGUn	0.2550	0.2918	0.0873	0.1403	99.94

Table 4. Table IV : Simulation results using real data.

	Houston HI		Samson HI		Jasper Ridge HI
	${\text{NRMSE}}_{𝒀}$	Time [s]	${\text{NRMSE}}_{𝒀}$	Time [s]	${\text{NRMSE}}_{𝒀}$	Time [s]
FCLS	0.2470	2.56	0.0545	1.38	0.2057	1.52
PLMM	0.0713	663.25	0.0239	103.84	0.0553	220.84
ELMM	0.0171	38.30	0.0119	14.76	0.0278	27.08
GLMM	0.0016	48.53	0.0006	46.69	0.0019	86.33
DeepGUn	0.2355	259.61	0.0862	121.88	0.1094	209.64

Equations

$${\mathbf{y}}_n=\mathbf{M}{\mathbf{a}}_n+{\mathbf{e}}_n,\text{subject to }{1}^{\top }{\mathbf{a}}_n=1\text{ and }{\mathbf{a}}_n\ge 0$$

$${\mathbf{y}}_n=f({\mathbf{M}}_0,{\mathbf{\theta }}_n){\mathbf{a}}_n+{\mathbf{e}}_n$$

$$\log p(X)\ge$$

$$\displaystyle-KL\big{(}q_{\phi}(Z|X)\|p(Z)\big{)}$$

$$\underset{{\mathcal{G}}_{\theta }}{\min }\underset{{\mathcal{C}}_{\varphi }}{\max }$$

$$\displaystyle+\operatorname{E}_{p(Z)}\big{\{}(1-\mathcal{C}_{\phi}(\mathcal{G}_{\theta}(Z)))\big{\}}\,.$$

$$\displaystyle\boldsymbol{m}_{p}\in\big{\{}\mathcal{G}_{\theta_{p}}(\boldsymbol{z}_{p})\,:\,\boldsymbol{z}_{p}\in\amsmathbb{R}^{K}\big{\}}$$

$${\mathbf{Z}}_0$$

$$\displaystyle{}={}\big{[}\mathcal{D}_{\phi_{1}}(\boldsymbol{m}_{1,0}),\ldots,\mathcal{D}_{\phi_{P}}(\boldsymbol{m}_{P,0})\big{]}\,.$$

$$J(\mathbf{A},\mathbb{Z})=\frac{1}{2}\sum _{n=1}^N\Vert {\mathbf{y}}_n-\tilde{\mathcal{G}}({\mathbf{Z}}_n){\mathbf{a}}_n{\Vert }_F^2+\mathcal{R}(\mathbf{A})+\mathcal{R}(\mathbb{Z})$$

$$\widetilde{\mathcal{G}}(\boldsymbol{Z}_{n})=\big{[}\mathcal{G}_{\theta_{1}}(\boldsymbol{z}_{1,n}),\ldots,\mathcal{G}_{\theta_{P}}(\boldsymbol{z}_{P,n})\big{]},\,\,\,n=1,\ldots,N$$

$$\mathcal{R}(\mathbf{A})=$$

$$\displaystyle\mathcal{S}^{1}{}={}\big{\{}\boldsymbol{A}\in\amsmathbb{R}^{P\times N}\,:\,\boldsymbol{A}\geq 0,\,\boldsymbol{1}^{\top}\boldsymbol{A}=\boldsymbol{1}^{\top}\big{\}}\,.$$

$$\mathcal{R}(\mathbb{Z})=\frac{{\lambda }_Z}{2}\sum _{n=1}^N\Vert {\mathbf{Z}}_n-{\mathbf{Z}}_0{\Vert }_F^2$$

$$(\hat{\mathbf{A}},\hat{\mathbb{Z}})={\mathrm{argmin}}_{\mathbf{A},\mathbb{Z}}J(\mathbf{A},\mathbb{Z}).$$

$$\begin{split}\min_{\mathbb{Z}}\,\,\,\frac{1}{2}\sum_{n=1}^{N}\Big{(}\|\boldsymbol{y}_{n}-\widetilde{\mathcal{G}}(\boldsymbol{Z}_{n})\boldsymbol{a}_{n}\|_{F}^{2}+\lambda_{Z}\|\boldsymbol{Z}_{n}-\boldsymbol{Z}_{0}\|_{F}^{2}\Big{)}\end{split}$$

$$\begin{array}{rl}\underset{\mathbf{A}}{\min }& \frac{1}{2}\sum _{n=1}^N\Vert {\mathbf{y}}_n-\tilde{\mathcal{G}}({\mathbf{Z}}_n){\mathbf{a}}_n{\Vert }_F^2+{\iota }_{{\mathcal{S}}^1}(\mathbf{A})\\ & +{\lambda }_A(\Vert {\mathcal{H}}_h(\mathbf{A}){\Vert }_{2,1}+\Vert {\mathcal{H}}_v(\mathbf{A}){\Vert }_{2,1}).\end{array}$$

$${\text{NRMSE}}_{\mathbb{X}}=\sqrt{\frac{\Vert \mathbb{X}-\hat{\mathbb{X}}{\Vert }_F^2}{\Vert \mathbb{X}{\Vert }_F^2}}.$$

$$\text{SAM}_{\mathbb{M}}=\frac{1}{N}\sum_{n=1}^{N}\sum_{p=1}^{P}\arccos\bigg{(}\frac{\boldsymbol{m}_{p,n}^{\top}\widehat{\boldsymbol{m}}_{p,n}}{\|\boldsymbol{m}_{p,n}\|\|\widehat{\boldsymbol{m}}_{p,n}\|}\bigg{)}.$$

Full text

I Introduction

Hyperspectral image analysis consists in a vast collection of algorithms and methods used to retrieve vital information from hyperspectral images (HI) in a increasing number of applications. Common applications include [1, 2] space exploration, remote sensing, surveillance, and, more recently, medical applications such as disease diagnosis and image-guided surgery [3]. One analysis methodology of particular interest is spectral unmixing (SU), which aims at retrieving sub-pixel information concerning the spectra of materials present in the scene, as well as estimating the proportions in which they contribute to each HI pixel [4, 5].

Many parametric models have been proposed to describe the interaction between light and the target surface [4, 6]. 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 pure material reflectances, 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 [6, 7, 8, 9] or variations of the EMs along the image [10, 11, 12]. More sophisticated models are required when such nonidealities have important impact on the formation of the HI.

I-A EM variability and learning-based SU methods

The variation of the endmembers across an HI (also called EM variability) is a very common effect since we can often associate multiple, different spectral signatures to each pure underlying material in a scene. EM variability can originate from environmental conditions, illumination, or atmospheric or temporal changes [13]. Its occurrence may incur the propagation of significant estimation errors throughout the unmixing process [10]. Different strategies have been proposed to cope with EM variability in SU. They can be classified in methods that represent EMs as sets, methods that model EMs as statistical distributions, and methods that incorporate parametric representations of EM variability in the mixing model [14].

Parametric models are raising considerable interest since they lead to good unmixing results and avoid the main drawbacks of the other groups of SU methods that address EM variability, namely the dependence on a priori knowledge of libraries of material spectra or the need for strong assumptions on the statistical distribution of the EMs for mathematical tractability [13, 14]. Recently proposed parametric models attempt to capture spectral variability by extending the LMM using either additive [10] or multiplicative [11, 12, 15, 16] scaling factors, or by considering tensor-based formulations [17, 18].

Although SU methods based on extended parametric models offer different trade-offs between representation capacity and model complexity, they still fail to achieve a desirable balance between a low-dimensional representation and enough flexibility to represent complex EM variability. Specifically, they fail to properly explore the fact that, although being very complex and spectrally non-homogeneous, spectral variability in real scenes is often confined to low-dimensional manifolds [19, 20, 21]. This property is due to the fact that the spectral signature of many materials is a function of only a few photometric or chemical properties of the medium. Prominent examples include packed particle spectra as a function of its roughness, size and density [19], leaf reflectance spectra as a function of various biophysical parameters [20], and soil reflectance as a function of moisture conditions [21]. Thus, existing models tend to be either too restrictive in their modeling capability or to lead to severely ill-posed estimation problems.

SU considering EM variability has also been formulated as a supervised learning problem, which is then solved without the need for an accurate physical model using neural networks (NNs) or support vector machines (SVMs) [22, 23, 24, 25]. However, these strategies depend on the availability of vast amounts of training data to adequately capture the spectral diversity of real scenes. This makes the training process computationally intensive and often intractable for large EM libraries, which must also be known a priori. Some works attempt to reduce the computational cost of these solutions by modifying learning algorithms to use hybrid soft-hard classification [26, 27, 28]. However, the resulting reconstructed abundance fractions do not have a clear physical interpretation due to the lack of a direct relationship to a physically motivated mixing model.

More recently, unsupervised SU approaches have also emerged by using autoencoders (AEC), which consist of encoder-decoder structured NNs originally devised for nonlinear dimensionality reduction [29]. These methods attempt to associate the decoder structure of the network with the LMM and the low-dimension representation of the input spectral vectors to the fractional abundances [30]. Different variations have been proposed, using pre-processing steps to reduce noise and outliers [31, 32], untiying the decoder from the encoder weights [33], using spectral angle distances to address nonlinear SU [34], or using denoising autoencoders to generate a robust initialization to matrix factorization-based SU strategies [35]. In [36] the authors proposed a nonlinear encoder-decoder structure to address the unmixing problem considering spectral variability. The proposed solution involves the simultaneous training of six neural networks to optimize a very large number of parameters. An autoencoder structure is employed to estimate the parameters of a spectral model of a hyperspectral image by minimizing the image reconstruction error while limiting the energy of some of the model parameters. No regularization strategy connecting the different pixels of the image is employed.

Despite their popularity, supervised learning-based SU algorithms are still not able to properly address the spectral variability problem, as they depend on extremely large amounts of labeled training data, leading to a computationally unfeasible learning process. Furthermore, the lack of a clear connection between AEC-based strategies and the physical mixing process makes one skeptic when concerning the robustness of AEC-based SU in face of more complex phenomena such as spectral variability.

I-B Proposed methodology

In this work, we propose a novel SU formulation that leverages the advantages of deep learning methods to address EM variability while still maintaining a strong connection to the physical mixing process, and using limited amounts of training data. Specifically, we adopt a deep generative NN to represent the manifold of EM spectra, which is then incorporated within the LMM. Generative models such as variational autoencoders (VAE) [37] and generative adversarial networks (GAN) [38] have recently obtained excellent performance at learning the probability distribution of complex data sets in very high dimensional spaces (e.g. natural images) from relatively small amounts of training data. The structure of generative models allow one to find a low-dimensional latent representation that parsimonously describes the variability of complex high-dimensional data sets. This leads to a low-dimensional parametrization of the training data distribution.

We formulate a novel unmixing strategy that can be cast as the problem of estimating the latent representations of the generative endmember models and the corresponding fractional abundances for each pixel in the HI. Specifically, we break down the SU problem in two steps. In the first step, we learn the latent EM variability manifold for each material in the scene using a deep generative EM model. The learning process uses pure pixel information directly extracted from the observed HI, which makes the proposed strategy suitable for unsupervised SU. In the second step, an alternating least-squares strategy is employed to estimate the parameters of an extended version of the LMM parametrized using the generative EM models obtained in the first step. The corresponding optimization problem is solved iteratively with respect to the abundances and to the low-dimensional representations of the EMs in the latent space of the deep generative models.

As a result, the proposed approach benefits from the reduced dimension of the latent space. Moreover, unlike current approaches, the new method does not depend on the careful selection of regularization parameters to yield a good performance. The resulting algorithm is named Deep Generative Unmixing algorithm (DeepGUn). The proposed method is strongly related to parametric models and leverages the learning and generalization capability of deep neural networks to properly represent the manifold of EM variability. Hence, DeepGUn leads to a model that is both low-dimensional and physically accurate, better describing the variability actually present in the scene.

Experimental results performed with both synthetic and real data indicate that the proposed strategy leads to more accurate abundance estimations than standard state-of-art SU methods accounting for EM variability. Qualitative analysis of the estimated abundance maps confirms these results. The improved accuracy comes at the expense of a small increase in the computational cost when compared to the best competing strategies.

This paper is organized as follows. Section II briefly reviews the LMM and its parametric extended versions. Section III discusses the basic properties of generative models in the context of VAE and GAN. Section IV introduces the proposed generative EM model and its learning strategy. In Section V we formulate the resulting SU problem, present the DeepGUn algorithm, and discuss aspects of the proposed optimization strategy. The neural network architecture is discussed in Section VI. The performance of the proposed method is compared with that of competing algorithms in Section VII. Finally, the conclusions are presented in Section VIII.

II Linear Mixing Models

The Linear Mixing Model (LMM) [4] assumes that a given $n$-th pixel $\boldsymbol{y}_{n}\in\amsmathbb{R}^{L}$, with $L$ bands, is represented as

[TABLE]

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

Different parametric models have been recently proposed to account for variable EM spectra within a given scene [11, 10, 12]. These models can be generically described as

[TABLE]

where $f$ is a parametric function, $\boldsymbol{M}_{0}\in\amsmathbb{R}^{L\times P}$ is a reference EM matrix, and $\boldsymbol{\theta}_{n}$ is a vector of parameters describing the manifold of EM variability.

Different functional forms have been proposed for $f(\boldsymbol{M}_{0},\boldsymbol{\theta}_{n})$ to account for EM variability in this framework, such as additive [10] or multiplicative [11, 12] variability factors acting upon the reference EM matrix $\boldsymbol{M}_{0}$. However, these models fail to achieve a desirable balance between a low-dimensional representation and enough flexibility to represent complex variability patterns. They tend to be either too restrictive in their modeling capability, or to lead to ill-posed optimization problems [18]. Instead of using a pre-defined parametric model, we propose to address this issue by learning a parametric function $f(\boldsymbol{M}_{0},\boldsymbol{\theta}_{n})$ using a generative model.

III Generative Models

Generative models attempt to estimate the probability distribution $p(X)$ of a random variable $X\in\amsmathbb{R}^{L}$ based on a set of observations $\boldsymbol{x}_{i}$, $i=1,\ldots,N_{x}$ in such a way that allows one to generate new samples that look similar to new realizations of $X$. The main characteristic of this problem, which sets it apart of other unsupervised learning methods such as density estimation, is the fact that we must to be able to sample from the estimated model $\hat{p}(X)$.

In many practical applications of interest, the dimensionality $L$ of the variable of interest $X$ is very high. This makes the general problem very difficult, as it amounts to estimating and sampling from an arbitrary high-dimensional probability density function [39, 40]. Nonetheless, the distributions of interest are often supported at a low-dimensional manifold of a set of so-called latent variables, and this fact can be explored to make the problem more tractable. A convenient way to address this problem is to define a new random variable $\amsmathbb{R}^{K}\ni Z\sim p(Z)$ with a known distribution in a low-dimensional space (e.g. an isotropic Gaussian distribution with $K\ll L$), and a parametric function (e.g. a neural network) $\mathcal{G}_{\theta}$ mapping $Z\mapsto\widehat{X}\in\amsmathbb{R}^{L}$ such that the image of $Z$ by $\mathcal{G}_{\theta}$ is a random variable whose distribution is very close to $p(X)$. In other words, the goal becomes to learn the parameters $\theta$ of $\mathcal{G}_{\theta}$ such that the distribution of $\widehat{X}=\mathcal{G}_{\theta}(Z)$ is as close to $p(X)$ as possible. Then, samples of $\widehat{X}$ can be generated by sampling from $Z\sim p(Z)$ and using the mapping $\mathcal{G}_{\theta}(Z)$.

Although estimating $\theta$ may still seem difficult at first, recent advances in machine learning such as VAEs [37] and GANs [38] have shown formidable performance at learning complex distributions such as those of natural images.

VAEs address this problem by assuming that the distribution of the observed data $X$ follows a directed graphical model $p(X|Z)$, which is represented by the function $\mathcal{G}_{\theta}$. The parameters of $\mathcal{G}_{\theta}$ are learned by maximizing a lower bound on the log-likelihood of $p(X)$ [37]:

[TABLE]

where $KL(\cdot\|\cdot)$ is the Kullback-Leibler divergence between two distributions, $\operatorname{E}_{\varsigma}\{\cdot\}$ is the expected value operator with respect to the distribution $\varsigma$ and $q_{\phi}(Z|X)$ is a variational approximation to the intractable posterior $p(Z|X)$, which is modeled by a function $\mathcal{D}_{\phi}$ (e.g. another neural network) parameterized in $\phi$. Note that $q_{\phi}(Z|X)$ must be a high-capacity distribution111Capacity of a distribution is a generic term to describe how complex a relationship it can model., so that it can provide a good approximation of the posterior $p(Z|X)$, which then allows the lower bound in (III) to be close to the true value of $\log\,p(X)$ [41].

GANs, on the other hand, attempt to learn the distribution $p(X)$ by searching for the Nash equilibrium of a two-player adversarial game [38]. A generator network $\mathcal{G}_{\theta}$ tries map the distribution of the latent variables $Z$ into the data distribution of $X$, and a discriminator network $\mathcal{C}_{\phi}$ tries to predict the probability of a random sample $\boldsymbol{x}_{i}$ coming from the true distribution $p(X)$ instead of being generated through $\mathcal{G}_{\theta}$. The generator $\mathcal{G}_{\theta}$ is trained by maximizing the probability of the discriminator making a mistake. This is formulated as the minimax optimization problem

[TABLE]

GANs are more flexible and have shown better performance at approximating complex distributions such as natural images (leading to sharper results) when compared to VAEs [38]. However, GANs are also much harder to train [40]. Moreover, VAEs naturally offer a way to obtain the latent representations corresponding to samples $\boldsymbol{x}_{i}\sim p(X)$ by mapping $X\mapsto Z$ using the function $\mathcal{D}_{\phi}$, which is also called an encoder model. This property and their stable training have motivated us to use VAEs in this work.

IV A Deep Generative Endmember Model

In this section, we propose to model the distribution of EM spectral variability using a deep generative model. By doing so, we can explicitly explore a common property of spectral variability: the EM spectra are usually confined to a low-dimensional manifold. This property is due to the fact that the spectral signature of many materials is a function of a few photometric or chemical properties of the medium. Prominent examples include packed particle spectra as a function of its roughness, size and density [19], leaf reflectance spectra as a function of various biophysical parameters [20], and soil reflectance as a function of moisture conditions [21].

IV-A The steps of the proposed SU method

We assume the existence of nonlinear functions $\mathcal{G}_{\theta_{p}}$, $p=1,\ldots,P$ (the generative model) that map latent representations $\boldsymbol{z}_{p}$ into their corresponding spectral signatures $\boldsymbol{m}_{p}$. We assume also the existence of encoder models $\mathcal{D}_{\phi_{p}}$ that map spectral signatures into their latent representations. In other words, we assume that any arbitrary observation $\boldsymbol{m}_{p}$ of a spectral signature of a material belongs to the set

[TABLE]

and thus can be equivalently represented by a corresponding low-dimensional vector $\boldsymbol{z}_{p}\in\amsmathbb{R}^{K}$ in the latent space of the generative model $\mathcal{G}_{\theta_{p}}$. This reasoning is illustrated in Fig. 1, where the encoder function $\mathcal{D}_{\phi}$ maps the input EM signature to the low-dimensional manifold. Reciprocally, low-dimensional vectors in the latent space can be mapped (decoded) to their corresponding spectral signatures using $\mathcal{G}_{\theta}$.

As such, we can formulate EM estimation in the SU problem in the latent domain (as opposed to the input spectral space), which is of a much lower order. Moreover, this approach will keep the physical interpretation of the model, provided that we have relevant training data to learn the generative models [42, 43, 44, 45]. This strategy relies on the existence of a priori training data for each material in the image, which might come in the form of, e.g., spectral libraries of laboratory measurements [46]. Nevertheless, we propose a more practical and effective approach to train the generators $\mathcal{G}_{\theta_{p}}$ by exploring information contained in multiple pure pixels extracted from the observed HI. The presence of multiple pure pixels in an observed HI is a characteristic of many real scenes, and can be leveraged to help in estimating the EM models, thus, reducing the ill-posedness of the SU problem222Pure pixels are 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..

Therefore, 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 methods, learn the generative and encoder models, $\mathcal{G}_{\theta_{p}}$ and $\mathcal{D}_{\phi_{p}}$, for all EMs in the scene ($p=1,\ldots,P$).

• ii)

Using the learned generative models, solve the SU problem by estimating the latent EM representations $\boldsymbol{Z}_{n}=[\boldsymbol{z}_{1,n},\ldots,\boldsymbol{z}_{P,n}]$ and the fractional abundance vectors $\boldsymbol{a}_{n}$ that can best represent the observed hyperspectral data, for all pixels in the scene ($n=1,\ldots,N$).

IV-B Learning the generative and encoder models $\mathcal{G}_{\theta_{p}}$ and $\mathcal{D}_{\phi_{p}}$

The objective of this first problem is to estimate the generative and encoding models $\mathcal{G}_{\theta_{p}}$ and $\mathcal{D}_{\phi_{p}}$, for $p=1,\ldots,P$. We assume the knowledge of a set $\mathcal{N}_{\mathcal{P},p}$ of pure pixels for the $p$-th EM, for all $p=1,\ldots,P$. Multiple pure pixels exist in many scenes, and can be directly extracted from the observed HI using automated EM extraction techniques [47, 48]. The sets of pure pixels $\mathcal{N}_{\mathcal{P},p}$, which can be seen as observations from the statistical distribution of each EM, are then used in the form of training data to learn the models $\mathcal{G}_{\theta_{p}}$ and $\mathcal{D}_{\phi_{p}}$ using a VAE [37]. If the set $\mathcal{N}_{\mathcal{P},p}$ is representative of the variability of the $p$-th material, the learned generative model $\mathcal{G}_{\theta_{p}}$ will be able to accurately describe the manifold of the $p$-th EM variability. Doing the same for all $p=1,\dots,P$ yields a set of variability models for all the EM spectra.

Although the extraction of multiple pure pixels from observed HIs is a well-established technique used to produce EM libraries [47], mixed pixels can sometimes be mistakenly identified as a pure pixel of some of the EMs. This constitutes a problem for library-based SU applications (e.g. MESMA and sparse SU) since some of the library spectra may end up not being representative of their EM class (material).

The smooth nature of the latent representation of VAEs allows the mitigation of this problem in the proposed approach. Assuming the availability of a reference EM matrix $\boldsymbol{M}_{0}$ of correctly identified signatures (which can be obtained using any EM extraction method) and of a set of encoder models $\mathcal{D}_{\phi_{p}}$, we can compute the latent representation of these reference signatures of each EM as

[TABLE]

where $\boldsymbol{m}_{p,0}$ is the $p$-th column of $\boldsymbol{M}_{0}$. The latent representation $\boldsymbol{z}_{p,0}$ can be used as a reference latent code for the $p$-th material. Thus, we can measure how close an estimated EM latent representation $\boldsymbol{z}_{p}$ is to the latent representation of a pure pixel by evaluating its Euclidean distance to $\boldsymbol{z}_{p,0}$. This can be performed since the output of VAEs have been shown to vary smoothly with changes of the latent variable [37]. Thus, we can use $\boldsymbol{Z}_{0}$ to regularize the SU problem to prevent $\mathcal{G}_{\theta_{p}}(\boldsymbol{z}_{p})$ from representing mixed pixels. This increases the robustness of the proposed approach.

IV-C Extracting sets of pure pixels from the observed HI

An important part of the proposed methodology consists in the extraction of the sets of pure pixels $\mathcal{N}_{\mathcal{P},p}$, $p=1,\ldots,P$ from the observed hyperspectral image $\boldsymbol{Y}$. Although different strategies have been proposed for image-based library construction (see e.g. [49, 50]), these techniques depend on multiple parameters that must be carefully adjusted in order to obtain good results. Instead of these approaches, we adopt a very simple strategy to select pure pixels from an HI that makes use of the reference matrix $\boldsymbol{M}_{0}$ extracted from the image using a pure-pixel-based endmember extraction algorithm (e.g. VCA [51]), which will also later be used to construct $\boldsymbol{Z}_{0}$ in (IV-B). We simply select as the elements of $\mathcal{N}_{\mathcal{P},p}$ the $S_{p}$ image pixels that have the smallest spectral angle to the reference signature in the $p$-th column of $\boldsymbol{M}_{0}$, where $S_{p}$ is the cardinality of $\mathcal{N}_{\mathcal{P},p}$ for $p=1,\ldots,P$. Although the success of this strategy depends on having a reasonably accurate estimation of $\boldsymbol{M}_{0}$, we experimentally found it to be more robust and easier to adjust than, for instance, the one in [47].

V The Unmixing Algorithm

Given a set of generative models $\mathcal{G}_{\theta_{p}}:\amsmathbb{R}^{K}\to\amsmathbb{R}^{L}$, $p=1,\ldots,P$ for each EM in the scene, a latent space representation $\boldsymbol{Z}_{0}$ of a reference EM matrix $\boldsymbol{M}_{0}$, and an HI $\boldsymbol{Y}=[\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{N}]$, the SU problem can be cast as the minimization of a risk functional of the form

[TABLE]

where $\boldsymbol{A}=[\boldsymbol{a}_{1},\dots,\boldsymbol{a}_{N}]\in\amsmathbb{R}^{P\times N}$ is the abundance matrix, $\mathbb{Z}\in\amsmathbb{R}^{N\times P\times K}$ is a 3-D tensor obtained by stacking all pixel-dependent latent EM representations $\boldsymbol{Z}_{n}$, such that $[\mathbb{Z}]_{n,:,:}=\boldsymbol{Z}_{n}$, $\mathcal{R}(\boldsymbol{A})$ and $\mathcal{R}(\mathbb{Z})$ are regularization terms to improve the problem conditioning, and the matrix-valued function $\widetilde{\mathcal{G}}(\boldsymbol{Z}_{n})$ defined as

[TABLE]

is the concatenation of the generative functions for each EM.

The term $\mathcal{R}(\boldsymbol{A})$ is a regularization functional that aims to provide spatial smoothness and to enforce positivity and sum-to-one constraints to the abundances. It is given by [12]

[TABLE]

where parameter $\lambda_{A}$ controls the contribution of this term to the cost function. The first two terms are a spatial regularizers over $\boldsymbol{A}$, where $\mathcal{H}_{h}$ and $\mathcal{H}_{v}$ are linear operators that compute the first-order horizontal and vertical gradients of a bidimensional signal, acting separately for each material of $\boldsymbol{A}$, and $\|\cdot\|_{2,1}$ is the $\mathcal{L}_{2,1}$ norm, defined as $\|\boldsymbol{X}\|_{2,1}=\sum_{n=1}^{N}\|\boldsymbol{x}_{n}\|_{2}$. The term $\iota_{\mathcal{S}^{1}}(\boldsymbol{A})$ is the indicator function of the unity simplex, i.e. $\iota_{\mathcal{S}^{1}}(\boldsymbol{A})=0$ if $\boldsymbol{A}\in\mathcal{S}^{1}$ and $\iota_{\mathcal{S}^{1}}(\boldsymbol{A})=\infty$ otherwise, where

[TABLE]

The term $\mathcal{R}(\mathbb{Z})$ constrains the EM latent representations $\mathbb{Z}$ to be close to the latent representation $\boldsymbol{Z}_{0}$ of the reference EM matrix $\boldsymbol{M}_{0}$. It is given by

[TABLE]

where parameter $\lambda_{Z}$ controls the contribution of this term to the cost function. This regularization makes the estimation problem more robust to the selection of the training data $\mathcal{N}_{\mathcal{P},p}$ by assuring the closeness of the estimated latent codes $\mathbb{Z}$ and the representations of pure pixels of each class. However, it relies indirectly on the reference EM signatures $\boldsymbol{M}_{0}$ (which are extracted from the observed HI with endmember extraction algorithms) being adequate representatives of their material classes in order to provide a good performance.

The optimization problem then becomes

[TABLE]

The problem defined in (11) is non-smooth and non-convex if solved simultaneously with respect to both variables $\boldsymbol{A}$, and $\mathbb{Z}$. However, an approximate solution can be found by minimizing (11) iteratively with respect to each variable, leading to the Deep Generative Unmixing (DeepGUn) method described in Algorithm 1. The DeepGUn algorithm consists of two distinctive steps. First, the generative endmember models generative and encoder models $\mathcal{G}_{\theta_{p}}$, $\mathcal{D}_{\phi_{p}}$, $p=1,\ldots,P$ are trained based on the pure pixels $\mathcal{N}_{\mathcal{P},p}$, $p=1,\ldots,P$ extracted from the observed HI and $\boldsymbol{Z}_{0}$ is computed. Afterwards, the alternating minimization approach is applied to compute the abundance maps and the latent representations of the endmembers for each pixel. We next describe the details of each optimization step. Implementation details are described in Sections VI and VII.

V-A Optimization with respect to $\mathbb{Z}$

Rewriting (11) considering only the terms in (7) that depend on $\mathbb{Z}$, the problem becomes

[TABLE]

This is a regularized nonlinear least squares problem, which can be solved individually for each pixel $\boldsymbol{y}_{n}$. Thus, (12) can be decomposed into $N$ non-convex, nonlinear optimization problems with dimensionality $K\times P$ by denoting each summand in (12) by $\widetilde{J}^{\,(n)}$, $n=1,\ldots,N$. We solve each of those problems $\widetilde{J}^{\,(n)}$ using a quasi-Newton algorithm, described in Algorithm 2, which provides an efficient solution for high-dimensional functions $\widetilde{\mathcal{G}}$ [52].

Although problem (12) is generally non-convex, recent research [53] has proven that, under suitable assumptions on the generator network $\widetilde{\mathcal{G}}$, the problem of recovering the latent variable $\boldsymbol{Z}_{n}$ does not have any stationary point (e.g. local minima or saddle points) outside a small neighborhood of the desired solution and its negative scalar multiple. This indicates the existence of a favorable global geometry of (12).

Note that $\widetilde{\mathcal{G}}$ is not necessarily differentiable with respect to the latent representations $\boldsymbol{Z}_{n}$, which can make the optimization problem more challenging. Nonetheless, quasi-Newton algorithms show excellent performance at non-smooth problems [54], where convergence is generally observed as long as the line search procedure does not return a point at which the objective function is non-differentiable. This allows quasi-Newton algorithms to be directly applied to obtain approximate solutions to non-smooth problems with good computational efficiency [54, 55].

V-B Optimization with respect to $\boldsymbol{A}$

Restating (11) considering only the terms in (7) that depend on $\boldsymbol{A}$ leads to

[TABLE]

Since the latent variables $\boldsymbol{Z}_{n}$ are fixed, (13) consists of a SU problem with a pixel-dependent EM matrix and an edge-preserving spatial regularization. Although this problem is not separable with respect to each pixel in the image, the Alternating Direction Method of the Multipliers (ADMM) can be used to obtain an efficient solution [56]. The solution of (13) using the ADMM is well described elsewhere (e.g. [11]) and will thus be omitted here for conciseness.

VI Neural Network Architecture

As discussed before, we used a VAE [37] to learn the generative and encoder models $\mathcal{G}_{\theta_{p}}$ and $\mathcal{D}_{\phi_{p}}$ from the sets of pure pixels $\mathcal{N}_{\mathcal{P},p}$. Compared to GANs, the training of VAEs is much simpler and more stable [40]. Moreover, VAEs naturally return the encoder model $\mathcal{D}_{\phi_{p}}$ as an approximation to the posterior distribution when learning $\mathcal{G}_{\theta_{p}}$. We have selected a dimension $K=2$ for the latent space, as it was experimentally verified to be sufficient to adequately capture the variability of each single material in a scene.

For the network architectures, we selected the number of layers and neurons according to the autoencoder implementation in [57, 29], with three hidden layers using ReLU activations (defined as $ReLU(x)=\max(x,0)$) in the hidden layers, which are described in more detail in Tables I and II.

We found that this configuration led to spectrally smooth generated signatures, and was effective at generalizing well with small training sample sizes. We trained the network for 50 epochs with the Adam optimizer [58] in TensorFlow, using a batch optimization with mini-batch size equal to one third of the total amount of training data for each EM.

VII Experimental Results

In this section, simulation results using both synthetic and real data illustrate the performance of the proposed method. We compare the proposed DeepGUn method with the fully constrained least squares (FCLS), the PLMM [10], the ELMM [11], and the GLMM [12]. In all experiments, the VCA algorithm [51] was used to extract the reference EM matrix $\boldsymbol{M}_{0}$ from the observed HI and to initialize the different methods. The abundance maps of all methods were initialized using the results obtained by the FCLS algorithm. The sets $\mathcal{N}_{\mathcal{P},p}$ of pure pixels were constructed by selecting the $100$ image pixels $\boldsymbol{y}_{n}$ with the smallest spectral angles relative to the reference EMs in $\boldsymbol{M}_{0}$. We ran the alternating optimization process in Algorithm 1 for at most $10$ iterations or until the relative change of $\boldsymbol{A}$ and $\mathbb{Z}$ was less than $10^{-3}$. The iterative procedure in Algorithm 2 was run until the relative change of $\boldsymbol{z}_{i}$ was less than $10^{-3}$. The performances were evaluated using the Normalized Root Means Squared Error (NRMSE) between the estimated abundance maps ($\text{NRMSE}_{\boldsymbol{A}}$), between the EM matrices ($\text{NRMSE}_{\mathbb{M}}$) and between the reconstructed images ($\text{NRMSE}_{\boldsymbol{Y}}$). The NRMSE between a true, generic tensor $\mathbb{X}$ and its estimate $\widehat{\mathbb{X}}$ is defined as

[TABLE]

Note that for the case of $\text{NRMSE}_{\boldsymbol{Y}}$, the reconstructed image $\widehat{\boldsymbol{Y}}$ is given by $[\widehat{\boldsymbol{Y}}]_{:,n}=[\widehat{\mathbb{M}}]_{:,:,n}[\,\widehat{\!\boldsymbol{A}}]_{:,n}$, $n=1,\ldots,N$.

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

[TABLE]

where $\boldsymbol{m}_{p,n}$ and $\widehat{\boldsymbol{m}}_{p,n}$ are the true and the estimated signatures of the $p$-th endmember in the $n$-th pixel, respectively.

VII-A Synthetic data

To quantitatively compare the different algorithms, four synthetic datasets were created, namely Data Cubes 1–4 (DC1–DC4), with 70$\times$70 pixels (DC1) and 50$\times$50 pixels (DC2, DC3 and DC4). These datasets were built using three (DC1, DC2 and DC3) and five (DC4) 224-band EMs extracted from the USGS Spectral Library [59] and spatially correlated abundance maps, as depicted in the first row of Fig. 2.

Spectral variability of the EMs was imposed using four different models. For the DC1 datacube, we adopted the variability model used in [10], consisting of pixelwise multiplicative spectral factors given by random piecewise-linear functions. For DC2, the variability model of [12] was used, consisting of band dependent scaling factors that varied smoothly in both the spatial and spectral dimensions. For DC3, we considered a simple model introduced in [49, Section IV-A-1] to emulate errors in atmospheric compensation as a function of the viewing geometry given the direct and diffuse light on the scene, and the solar path transmittance. For datacube DC4, we used as endmembers pure pixels of five materials (asphalt, tree, roof, metal and dirt) which were manually extracted from a real hyperspectral image, thus depicting realistic spectral variability. White Gaussian noise was finally added to all datasets to yield a 30dB SNR.

The optimal parameters for each algorithm were selected by performing grid searches for each dataset. The ranges in which the parameters were searched were selected according to those discussed by the authors in the original publications. For the PLMM we searched for $\alpha$, $\beta$ and $\gamma$ in the ranges $[0.01,\,0.1,\,0.35,\,0.7,\,1.4,\,25]$, $[10^{-9},\,10^{-5},\,10^{-4},\,10^{-3}]$ and $[10^{-2},\,0.1,\,1,\,10,\,10^{2}]$, respectively. For both ELMM and GLMM, the parameters were selected among the following values: $\lambda_{S},\,\lambda_{M}\in[0.01,\,0.1,\,1,\,5,\,10,\,50]$, $\lambda_{A}\in[0.001,\,0.01,\,0.05]$, and $\lambda_{\psi},\,\lambda_{\mathbb{\Psi}}\in[10^{-6},\,10^{-3},\,1,\,10^{3}]$. For the proposed DeepGUn algorithm, we fixed $\lambda_{Z}=0.1$ and selected $\lambda_{A}$ among the values $[0.005,\,0.01,\,0.05]$. For the proposed method, the sets of pure pixels for each EM $\mathcal{N}_{\mathcal{P},p}$ were constructed by selecting the $100$ pixels closest to the reference materials $\boldsymbol{M}_{0}$.

The quantitative results are shown in Table III, with the best results for each metric marked in bold. The proposed method clearly outperformed the competing algorithms in terms of $\text{NRMSE}_{\boldsymbol{A}}$ for all four datasets. Qualitatively, the abundance maps provided by DeepGUn, displayed in Fig. 2, are clearly much closer to the true abundance maps than those provided by the other methods. These are important results, as accuracy in abundance estimation is the main objective of SU.

For the EM reconstruction metrics $\text{NRMSE}_{\mathbb{M}}$ and $\text{SAM}_{\mathbb{M}}$, DeepGUn gave the best results for all data cubes. This indicates that the proposed endmember model used by DeepGUn allows for precise material identification from the observed hyperspectral scenes.

The reconstruction error $\text{NRMSE}_{\boldsymbol{Y}}$ of the DeepGUn algorithm was comparable to the FCLS and significantly larger than that of the GLMM. This is natural since the GLMM has more degrees of freedom. However, the connection between $\text{NRMSE}_{\boldsymbol{Y}}$ and the abundance reconstruction error is far from being direct, as can be seen in Table III.

The execution times, at the rightmost column of Table III, indicate that the computational complexity of DeepGUn is somewhere between the complexities of GLMM and PLMM, the two major competing algorithms. Hence, the DeepGUn method yielded superior SU performance, with easier parameter tuning, and at a reasonable computational cost.

VII-B Influence of the latent dimension $K$

An important parameter in the design of the DeepGUn method is the dimensionality $K$ of the latent space of the generative endmember models $\mathcal{G}_{i}$, $i=1,\ldots,P$. The individual dimensions of the latent space are used to represent changes in the endmember signatures due to spectral variability. Since in a given scene the endmembers are likely to be affected only by a small number of effects (hence the hypothesis that they are supported at a low-dimensional manifold), $K$ should be small in order to avoid introducing spurious effects and increasing the computational complexity of the solution.

To illustrate this, we performed a simulation with DeepGUn where we varied the dimensionality of the latent space $K$ and measured the normalized abundance reconstruction error $\text{NRMSE}_{\boldsymbol{A}}$. For this, we considered the data cubes DC1 and DC2 described above in Section VII-A. The results are depicted in Figure 3, and show that the abundance estimation error tends to increase with $K$. It can also be seen that there is a sharper increase in $\text{NRMSE}_{\boldsymbol{A}}$ for DC1 when compared to DC2. This is likely due to the fact that the dataset DC2 uses a more complex model to generate endmember variability. This indicates that the selection of a small value for $K$ is important to obtain good unmixing results.

VII-C Real data

We considered the Houston, Samson and the Jasper Ridge datasets for the simulations with real data. These datasets were captured by the AVIRIS instrument, and originally had 224 bands. The spectral bands corresponding to water absorption and low SNR regions were removed, resulting in 188 bands for the Houston image, 156 bands for the Samson image and 198 bands for the Jasper Ridge image. Previous studies indicate that the Houston HI has four predominant EMs [11], while the Samson and Jasper Ridge HIs are known to have three and four EMs, respectively [60].

The reconstructed abundance maps for both datasets and all algorithms are shown in Figs. 4, 5 and 6. For the Houston dataset, the last row of Fig. 4 shows that the abundance maps provided by the DeepGUn method better evidence the strong vegetation and concrete abundances at the stadium field and stands, respectively, as well as the stronger asphalt abundances in the parking lot. For the Samson and Jasper Ridge images, a clear performance improvement can be seen for the DeepGUn algorithm. Note, for instance, a smaller confusion between the Water and Soil EMs in the Samson HI when compared to the other methods. Similarly, for the Jasper Ridge HI, the DeepGUn method leads to considerably stronger Water abundances in the region containing the river. Moreover, although the ELMM provided a better estimation of the road in the scene when compared to the remaining methods, it also resulted in a greater confusion between the Vegetation and Soil EMs, especially in the right part of the scene.

The quantitative results for all algorithms and datasets are shown in Table IV. Since the correct abundance values (the ground truth) are not available for most real images, the reconstruction error $\text{NRMSE}_{\boldsymbol{Y}}$ has been used as a sort of quality verification. As was the case for the synthetic data, the DeepGUn reconstruction errors are higher than those yielded by other methods that address spectral variability. However, the reconstruction error is definitely not a good performance measure for abundance estimation in real images, which is the main objective of unmixing algorithms. The higher reconstruction errors of DeepGUn in this case are just due the fact that DeepGUn has much fewer degrees of freedom than the ELMM, PLMM and GLMM algorithms. In fact, the DeepGUn has only $K\times P$ degrees of freedom for each pixel, which is comparable to the FCLS ($P$) much smaller than the ELMM, GLMM and PLMM methods ($>L\times P$). Although this means that the ELMM, GLMM and PLMM can achieve arbitrarily small reconstruction errors $\text{NRMSE}_{\boldsymbol{Y}}$, this is not necessarily reflected as good abundance estimation results. The execution times of the proposed DeepGUn method were again comparable to those of the other algorithms addressing spectral variability, which indicates that it scales well with larger image sizes.

VIII Conclusions

In this paper, a deep generative EM model was proposed to address spectral variability in SU of HIs. Instead of relying on user-defined parametric EM models which have shown to be very hard to estimate in practical scenes, the proposed methodology leveraged the generalization capability of deep neural networks to accurately model EM spectra while still maintaining a strong connection to the physical mixing process. A deep generative model for each EM was trained prior to unmixing by using pure pixel information extracted directly from the observed HI, which allowed for an unsupervised formulation. The proposed EM model was then applied to solve the SU problem, which was cast as the estimation of the low-dimensional representations of the EMs in the latent space of the deep generative models and their corresponding fractional abundances, for each pixel. The resulting DeepGUn algorithm presented excellent performance despite the simple strategy used for selecting the training data for learning the generative model. Simulations using synthetic and real data indicate that the proposed method can lead to significant improvements in abundance estimation accuracy.