Member, IEEE, Tales Imbiriba, Pau Closas, Senior Member, IEEE Ricardo A Borsoi
email: ricardo.borsoi@univ-lorraine.fr

Dynamical Hyperspectral Unmixing with Variational Recurrent Neural Networks

Keywords: Hyperspectral data, hyperspectral unmixing, recurrent neural networks, deep learning, multitemporal

Multitemporal hyperspectral unmixing (MTHU) is a fundamental tool in the analysis of hyperspectral image sequences. It reveals the dynamical evolution of the materials (endmembers) and of their proportions (abundances) in a given scene. However, adequately accounting for the spatial and temporal variability of the endmembers in MTHU is challenging, and has not been fully addressed so far in unsupervised frameworks. In this work, we propose an unsupervised MTHU algorithm based on variational recurrent neural networks. First, a stochastic model is proposed to represent both the dynamical evolution of the endmembers and their abundances, as well as the mixing process. Moreover, a new model based on a low-dimensional parametrization is used to represent spatial and temporal endmember variability, significantly reducing the amount of variables to be estimated. We propose to formulate MTHU as a Bayesian inference problem. However, the solution to this problem does not have an analytical solution due to the nonlinearity and non-Gaussianity of the model. Thus, we propose a solution based on deep variational inference, in which the posterior distribution of the estimated abundances and endmembers is represented by using a combination of recurrent neural networks and a physically motivated model. The parameters of the model are learned using stochastic backpropagation. Experimental results show that the proposed method outperforms state of the art MTHU algorithms.

I. INTRODUCTION

Hyperspectral images (HIs) have very high spectral resolution, which allows for a precise discrimination of different materials in a scene [START_REF] Bioucas-Dias | Hyperspectral remote sensing data analysis and future challenges[END_REF]. However, physical limitations of spectral image acquisition and large distances between the sensor and the scene of interest as seen in, e.g., remote sensing, means that each pixel of an HI may cover a large area of the scene and typically contains a mixture of different materials [START_REF] Shaw | Spectral imaging for remote sensing[END_REF]. Hyperspectral unmixing (HU) aims to decompose an HI into the spectral signatures of the pure materials it contains (the endmembers -EMs), and the proportions with which they appear in each pixel (the abundances) [START_REF] Keshava | Spectral unmixing[END_REF].

The classical approach to describe the interactions between light and the different materials in a pixel is the linear mixing model (LMM) [START_REF] Keshava | Spectral unmixing[END_REF]. However, the LMM assumes the EM signatures to be the same for all pixels in an HI, disregarding spectral variability of the EMs which can be caused by, e.g., atmospheric, illumination or seasonal variations, and propagates significant errors throughout the HU processing chain [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF], [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF]. Thus, significant effort has been dedicated to address spectral variability in HU (see Section II for a review).

More recently, multitemporal HU (MTHU) has been receiving increasing interest in the literature since it leverages information in sequences of HIs acquired at different time instants to reveal the dynamical evolution of the endmembers and abundances in a scene [START_REF] Henrot | Dynamical spectral unmixing of multitemporal hyperspectral images[END_REF]- [START_REF] Liu | Bayesian unmixing of hyperspectral image sequence with composite priors for abundance and endmember variability[END_REF]. MTHU has proven important for many applications such as invasive species mapping in rainforests [START_REF] Somers | Invasive species mapping in hawaiian rainforests using multi-temporal hyperion spaceborne imaging spectroscopy[END_REF], [START_REF]Multi-temporal hyperspectral mixture analysis and feature selection invasive species mapping in rainforests[END_REF], and monitoring vegetation cover in shrublands [START_REF] Lippitt | Multidate MESMA for monitoring vegetation growth forms in southern california shrublands[END_REF] or seasonal variations of vegetation cover in dry forests [START_REF] Goenaga | Unmixing analysis of a time series of hyperion images over the guánica dry forest in puerto rico[END_REF]. Moreover, MTHU is also useful to perform change detection at the subpixel level [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF], [START_REF] Guo | Change detection for hyperspectral images via convolutional sparse analysis and temporal spectral unmixing[END_REF]. However, spectral variability can be very significant in MTHU due to different seasonal and acquisition conditions [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF], [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF], [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF].

Addressing both the spatial and temporal spectral variability of the EMs is challenging, and has only been done in MTHU by supervised techniques [START_REF] Somers | Invasive species mapping in hawaiian rainforests using multi-temporal hyperion spaceborne imaging spectroscopy[END_REF], [START_REF] Borsoi | Fast unmixing and change detection in multitemporal hyperspectral data[END_REF]. However, supervised MTHU techniques require prior knowledge of libraries containing spectral signatures which can accurately represent the endmembers for each image in the time sequence. Such libraries can be difficult or expensive to collect. Unsupervised MTHU methods, on the other hand, estimate both the endmember signatures and the abundances for all time instants directly from the observed HI sequence. Thus, unsupervised methods are of great practical interest, but can be challenging to design. See Section II for a review of MTHU methods.

Machine learning has become a popular framework to solve the HU problem [START_REF] Bhatt | Deep learning in hyperspectral unmixing: A review[END_REF]. Recent developments include methods based on, e.g., autoencoders (AECs) [START_REF] Palsson | Blind hyperspectral unmixing using autoencoders: A critical comparison[END_REF] or unrolled optimization-based neural networks [START_REF] Zhou | ADMM-Based hyperspectral unmixing networks for abundance and endmember estimation[END_REF] (see Section II for a review). In particular, solutions based on deep learning are especially attractive for HU when the mixing model considers nonidealities such as, e.g., nonlinearity [START_REF] Wang | Nonlinear unmixing of hyperspectral data via deep autoencoder networks[END_REF] and EM variability [START_REF] Borsoi | Deep generative endmember modeling: An application to unsupervised spectral unmixing[END_REF], circumventing the need to construct complex analytical models to represent such physical effects.

However, the literature lacks MTHU solutions that are unsupervised and take spatial and temporal EM variability into account, which are addressed in this work. In particular, we also address several other needs, including the development of parsimonious models for EM variability with adjustable flexibility, and of machine learning-based strategies for MTHU which jointly leverage both a physically motivated and datadriven (e.g., neural networks) models in a principled manner. Such hybrid approaches, where physics-informed models are used to regularize and provide interpretability to data-driven methodologies are becoming increasingly popular [START_REF] Imbiriba | Hybrid Neural Network Augmented Physics-based Models for Nonlinear Filtering[END_REF].

In this work, we propose an unsupervised MTHU algorithm based on variational recurrent neural networks (RNNs). First, a stochastic model is proposed to represent both the dynamical evolution of the EMs and of the abundances, as well as the mixing process. Moreover, a new low-dimensional model is used to represent spatial and temporal EM variability by parametrizing band-dependent scaling variations of the EMs using a set of smooth spectral basis functions. This allows us to control the flexibility of the model by varying the number and types of basis functions. To model the abundances, we make use of the softmax basis representation [START_REF] Barber | Bayesian Reasoning and Machine Learning[END_REF], which leads to a physically accurate model and has been successfully used for fuzzy classification [START_REF] Kent | Spatial classification using fuzzy membership models[END_REF] and single-image HU [START_REF] Eches | Enhancing hyperspectral image unmixing with spatial correlations[END_REF]. This way, we can use a Gaussian distribution to represent the abundances in the softmax basis, which closely approximates a Dirichlet distribution when mapped back to the original abundance domain (i.e., the unit simplex) [START_REF] Mackay | Choice of basis for Laplace approximation[END_REF].

MTHU is then formulated as a Bayesian inference problem. However, exact inference is analytically intractable due to the nonlinearities in the model. Thus, we consider a variational inference solution based on RNNs, in which the approximate joint posterior distribution of abundances and EMs is learned by maximizing a lower bound over the marginal likelihood of all pixels. Note that approximating the true posterior typically requires a flexible family of distributions which can be represented using neural networks [START_REF] Kingma | Auto-encoding variational bayes[END_REF]. However, using feedforward neural networks leads to models with large numbers of parameters, making inference costly. By exploiting the temporal structure in the data (e.g., Markovity), RNNs provide a solution that gives flexibility while also having a lower number of parameters (being computationally lightweight). Besides, RNNs have shown excellent performance in numerous sequence modeling tasks [START_REF] Lipton | A critical review of recurrent neural networks for sequence learning[END_REF]. Interpretability of the estimated abundances and EMs is paramount for the applicability of MTHU systems. For this reason, we parameterized the joint posterior distribution using a hybrid model composed of physics-based and data-driven components. More specifically, the posterior is modeled by a family of nonlinear functions constructed by integrating both a simple, physically motivated model that is able to provide an approximate abundance estimate, and a bidirectional RNN that can represent more complex effects (i.e. not captured by the simpler model). The parameters of the model and of the posterior distribution are learned based on all image pixels using stochastic gradient descent (SGD). The contributions of this paper include:

' a new low-dimensional model to represent the spatial and temporal variability of the EMs with a small amount of parameters to be learned; ' a stochastic model describing the temporal evolution of the abundances and of the EM variability parameters, which leverages the softmax basis used in single-image HU [START_REF] Eches | Enhancing hyperspectral image unmixing with spatial correlations[END_REF] to obtain a physically accurate representation of the abundance dynamics using a Gaussian distribution; ' a deep variational inference

II. BACKGROUND AND RELATED WORK

A general multitemporal linear mixing model represents the n-th pixel of an HI acquired at time t as:

y n,t " M n,t a n,t `rn,t , s.t. 1 J a n,t " 1, a n,t ě 0 , (1)
where, for each time t P t1, . . . , T u and pixel n P t1, . . . , N u, y n,t P R L denotes the observed pixel with L bands, the columns of M n,t P R LˆP contain the spectral signatures of the P endmembers in the scene, vector a n,t P R P contains the fractional abundances of each EM, and r n,t P R L represents additive noise. Note that the general model ( 1) can accommodate EM variability both in space and in time, being able to represent effects such as, e.g., atmospheric, illumination, or seasonal variations [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF]. In the following, we review different HU strategies addressing multitemporal sequences, spatial EM variability, and based on deep learning frameworks.

A. Multitemporal HU

A fundamental aspect of MTHU is taking into account the relationship between the EMs and abundances at different time instants. Since these are usually temporally correlated, this can greatly improve the performance of unmixing algorithms. Most previous works have been focused on addressing the variability of the EMs in time. This was usually performed by considering a more constrained version of model [START_REF] Bioucas-Dias | Hyperspectral remote sensing data analysis and future challenges[END_REF], where only the temporal variability of the EMs is considered, leading to M n,t " M m,t , for all 1 ď n, m ď N [6]- [START_REF] Borsoi | Kalman filtering and expectation maximization for multitemporal spectral unmixing[END_REF].

Several works have considered parametric models to represent the temporal variability of the EMs (i.e., using only a single EM matrix per image). For instance, dynamical model was used in [START_REF] Henrot | Dynamical spectral unmixing of multitemporal hyperspectral images[END_REF] to constrained the EMs to be a scaled versions of a reference EM matrix, with smoothly varying scaling factors. Another model considered the EM matrices at each time instant to be an additive perturbation of a mean EM matrix [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF]. Using this model, MTHU was performed using a two-stage stochastic programming approach [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF]. Other works considering this model have also proposed MTHU solutions using a distributed algorithm using sparsity constraints [START_REF] Sigurdsson | Sparse distributed multitemporal hyperspectral unmixing[END_REF], and a hierarchical Bayesian framework incorporating additive residual terms to account for outliers [START_REF] Thouvenin | A hierarchical bayesian model accounting for endmember variability and abrupt spectral changes to unmix multitemporal hyperspectral images[END_REF]. A recent work proposed a hierarchical Bayesian MTHU strategy (called HBUN) which incorporated priors promoting the spectral smoothness of the estimated EM signatures, and the spatial and temporal smoothness of the abundances [START_REF] Liu | Bayesian unmixing of hyperspectral image sequence with composite priors for abundance and endmember variability[END_REF]. Another model representing the EMs using bandwise multiplicative scalings of a set of reference EM signatures was considered in [START_REF] Borsoi | Kalman filtering and expectation maximization for multitemporal spectral unmixing[END_REF]. MTHU was performed by combining a Bayesian filtering (the Kalman filter and smoother) with the expectation maximization method. Heteroscedastic measurement noise was also considered in [START_REF] Zhuo | A novel spectral-temporal Bayesian unmixing algorithm with spatial prior for Sentinel-2 time series[END_REF], where the (diagonal) covariance matrix of

KLD MSE Transition

Variational Posterior

Fig. 1: Illustrative diagram of the proposed ReDSUNN method. A likelihood and transition PDFs, with parameters θ, define the spatiotemporal mixing process (i.e., the generative model). The variational posterior, with parameters φ, approximates the unmixing solution (i.e., the PDF of the abundances and EMs conditioned on the HIs) recursively, being implemented using an RNN. The parameters of these PDFs are learned by maximizing a loss function based on the ELBO, which balances data reconstruction (MSE) and consistency between posterior and prior (KLD). The parameters cn,t and ψ n,t are mapped to the abundances and EM matrices, an,t and M n,t using a deterministic model. The notation z ´1 represents a delay. See Section III for more details. Fig. 2: In the generative model (left), conditional independence among pixels y t is assumed given abundances at and variability coefficients ψ t , as well as Markovity of state variables. During inference (right), state estimation exploits correlations between pixels, states and hidden representations ht generated by a bi-directional RNN.

the measurement noise was estimated along with the EMs and abundances in a maximum a posteriori framework.

However, these methods disregard EM variability within a single image, which limits their performance. Existing MTHU algorithms which tackle both spatial and temporal variabilty are based on the MESMA framework [START_REF] Roberts | Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[END_REF], which searches for a combination of EM signatures in a library which can best reconstruct a given HI, rendering it highly interpretable. These have been applied for MTHU in vegetation monitoring applications [START_REF] Somers | Invasive species mapping in hawaiian rainforests using multi-temporal hyperion spaceborne imaging spectroscopy[END_REF], [START_REF] Lippitt | Multidate MESMA for monitoring vegetation growth forms in southern california shrublands[END_REF], [START_REF] Dudley | A multi-temporal spectral library approach for mapping vegetation species across spatial and temporal phenological gradients[END_REF]. Recent advances were also made on developing efficient solutions to this problem with theoretical guarantees by exploring the temporal information of the abundances [START_REF] Borsoi | Fast unmixing and change detection in multitemporal hyperspectral data[END_REF]. However, MESMA-based approaches are computationally costly and their performance depends strongly on the accuracy of the library, rendering this model inadequate for unsupervised processing. Note that some works proposed to leverage complementary high-resolution (Landsat) images to unmix low-resolution (MODIS) images [START_REF] Wang | Spatio-temporal spectral unmixing of time-series images[END_REF], [START_REF]Real-time spatiotemporal spectral unmixing of MODIS images[END_REF]. However, the availability of such complementary information and difficulties associated to differences in their acquisition times limit the applicability of such methods in practice.

B. HU with EM variability

The variability of the EMs across an HI occurs due to atmospheric, illumination (e.g., topographic) or intrinsic variations of the EM spectra [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF], [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF], [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF]. It introduces errors which propagate through the various steps of the traditional HU processing chains and can have significant negative impact on the abundance estimation performance. EM variability is generally addressed in HU by representing the spectral signature of each material using structured spectral libraries, statistical distributions, or physically motivated parametric models [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF].

Spectral library-based methods represent the EMs in each pixel as one of several spectral signatures in a dictionary known a priori, which implies formulating HU as a structured sparse regression problem. The sparsity prior can be addressed either with combinatorial approaches, which are computationally costly [START_REF] Roberts | Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[END_REF], or using different relaxations of the L 0 seminorm based on convex [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF], [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF] or nonconvex [START_REF] Drumetz | Hyperspectral image unmixing with endmember bundles and group sparsity inducing mixed norms[END_REF] sparsity promoting penalties. Other strategies also allow each EM to be a convex combination of library spectra [START_REF] Uezato | Hyperspectral unmixing with spectral variability using adaptive bundles and double sparsity[END_REF]. Such relaxations are computationally easier to solve. This can be a reasonable modeling assumption when multiple signatures of the same EM can contribute to form a single pixel. However, such relaxations might reduce the interpretability of the solutions if the goal is to identify which signature is active [START_REF] Drumetz | Hyperspectral image unmixing with endmember bundles and group sparsity inducing mixed norms[END_REF]. Nevertheless, the performance of both kinds of strategies is strongly dependent on the quality of the spectral library.

Using statistical distributions to model the endmembers has been well-investigated as it provides principled HU solutions through a Bayesian framework. The Gaussian [START_REF] Halimi | Unsupervised unmixing of hyperspectral images accounting for endmember variability[END_REF], mixture of Gaussians [START_REF] Zhou | A gaussian mixture model representation of endmember variability in hyperspectral unmixing[END_REF] or Beta [START_REF] Du | Spatial and spectral unmixing using the beta compositional model[END_REF] distributions have been considered. However, when complex distributions (such as the Beta or mixtures of Gaussians) are used to represent the endmembers, HU (which consists in a Bayesian inference problem) can become computationally expensive.

Another approach to address spectral variability consists in representing the signatures of the EMs in each pixel using a physically meaningful parametric model, and estimating the model parameters during HU. Examples of such models include the use of additive perturbations [START_REF] Thouvenin | Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[END_REF], spectrally uniform or spectrally localized multiplicative scaling factors [START_REF] Drumetz | Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability[END_REF], [START_REF] Imbiriba | Generalized linear mixing model accounting for endmember variability[END_REF], or a combination thereof [START_REF] Hong | An augmented linear mixing model to address spectral variability for hyperspectral unmixing[END_REF]. Other models explicitly exploit spatial information, using multiscale [START_REF] Borsoi | A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF] or low-rank tensor representation [START_REF] Imbiriba | Low-rank tensor modeling for hyperspectral unmixing accounting for spectral variability[END_REF], and external information (e.g., LiDAR) by means of digital surface models [START_REF] Uezato | Illumination invariant hyperspectral image unmixing based on a digital surface model[END_REF]. However, designing physically accurate models whose parameters can at the same time be properly recovered from an HI can be challenging.

C. Deep learning-based HU

Deep learning has recently become a popular approach to perform HU. While HU was traditionally viewed as a regression problem (i.e., learning a mapping from the pixels to the abundances) [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF], [START_REF] Plaza | Spectral mixture analysis of hyperspectral scenes using intelligently selected training samples[END_REF], recent work has been focused on developing unsupervised or self-supervised strategies, which avoid the need for vast amounts of training data. Among such strategies, AECs have become a predominant approach for deep learning-based HU due to their close connection to linear or nonlinear mixing models and good experimental performance [START_REF] Palsson | Blind hyperspectral unmixing using autoencoders: A critical comparison[END_REF]. The latent representations of the image pixels obtained by the AEC are associated to the abundances, and the decoder network to the mixing model [START_REF] Guo | Hyperspectral image unmixing using autoencoder cascade[END_REF], [START_REF] Palsson | Hyperspectral unmixing using a neural network autoencoder[END_REF]. Several AEC methods for HU have been proposed for linear HU by using different choices of encoder networks including, e.g., denoising layers [START_REF] Su | DAEN: Deep autoencoder networks for hyperspectral unmixing[END_REF], [START_REF] Sahoo | Deep-learning-based latent space encoding for spectral unmixing of geological materials[END_REF], spatial-spectral (convolutional) architectures [START_REF] Palsson | Convolutional autoencoder for spectral-spatial hyperspectral unmixing[END_REF], [START_REF] Hua | Gated autoencoder network for spectral-spatial hyperspectral unmixing[END_REF] and the use of sparsity constraints [START_REF] Qu | uDAS: An untied denoising autoencoder with sparsity for spectral unmixing[END_REF].

AECs have also been used to perform nonlinear HU by considering nonlinear decoder networks to address complex mixing effects. This includes a post-nonlinear mixing model [START_REF] Wang | Nonlinear unmixing of hyperspectral data via deep autoencoder networks[END_REF], additive nonlinearities [START_REF] Zhao | Hyperspectral unmixing for additive nonlinear models with a 3-D-CNN autoencoder network[END_REF] and the use of application-specific nonlinear neural network layers [START_REF] Shahid | Unsupervised hyperspectral unmixing via nonlinear autoencoders[END_REF]. Another work also exploited the relationship between the encoder and decoder networks to propose a model-based architecture [START_REF] Li | Model-based deep autoencoder networks for nonlinear hyperspectral unmixing[END_REF].

Spectral variability was also addressed using deep learning methods. In [START_REF] Borsoi | Deep generative endmember modeling: An application to unsupervised spectral unmixing[END_REF], a generative EM model is proposed to represent the variability of the EMs on a low-dimensional manifold. Such a model was used to perform HU using methods inspired by matrix factorization [START_REF] Borsoi | Deep generative endmember modeling: An application to unsupervised spectral unmixing[END_REF], sparse regression [START_REF] Borsoi | Deep generative models for library augmentation in multiple endmember spectral mixture analysis[END_REF] and probabilistic approaches [START_REF] Sahoo | Deep-learning-based latent space encoding for spectral unmixing of geological materials[END_REF], [START_REF] Shi | Probabilistic generative model for hyperspectral unmixing accounting for endmember variability[END_REF]. Gaussian process regression has also been considered as a non-parametric approach to mitigate the effects of spectral variability in HU [START_REF] Uezato | A novel spectral unmixing method incorporating spectral variability within endmember classes[END_REF], [START_REF] Koirala | Robust supervised method for nonlinear spectral unmixing accounting for endmember variability[END_REF].

Other approaches also used multiple sets of pure pixels extracted in a self-supervised manner to regularize (explicitly or implicitly) AEC-based HU algorithms in order to improve their robustness [START_REF] Jin | Tanet: An unsupervised two-stream autoencoder network for hyperspectral unmixing[END_REF]- [START_REF] Li | Self-supervised robust deep matrix factorization for hyperspectral unmixing[END_REF]. An approach to learn learn dynamical models for the spectra of individual materials has also been proposed [START_REF] Drumetz | Learning endmember dynamics in multitemporal hyperspectral data using a statespace model formulation[END_REF]. Finally, other methods have been proposed using, e.g., Wasserstein [START_REF] Min | JMnet: Joint metric neural network for hyperspectral unmixing[END_REF], adversarial [START_REF] Jin | Adversarial autoencoder network for hyperspectral unmixing[END_REF], or cycle-consistency [START_REF] Gao | CyCU-Net: Cycle-consistency unmixing network by learning cascaded autoencoders[END_REF] loss functions during training. Deep priors [START_REF] Rasti | UnDIP: Hyperspectral unmixing using deep image prior[END_REF] and unfolding optimization-based neural networks have also been considered [START_REF] Zhou | ADMM-Based hyperspectral unmixing networks for abundance and endmember estimation[END_REF].

We highlight, however, that despite these advances previous MTHU methods did not address spatial EM variability without supervision, and also did not exploit deep learning frameworks. In the following, we will present a probabilistic model representing both spatial and temporal EM variability, and develop an unsupervised inference strategy based on RNNs.

III. OVERVIEW OF THE PROPOSED APPROACH

We consider a probabilistic framework for MTHU. This amounts to two steps. The first is the modeling step, which consists in defining a set of probability density functions (PDFs) which describe how the abundances and EMs evolve over time, and how the pixels are generated. Note that the model PDFs typically depend on different deterministic parameters, some of which are specified a priori (which we refer to as hyperparameters) and some which we intend to learn from the observed HIs. We represent the parameters to be learned in the set θ. We include the subscript θ on the model PDFs in order to make their dependence on these parameters explicit.

The second is the inference step, which consists in computing (an approximation of) the posterior distribution, which is the PDF of the abundances and EMs conditioned on the observed pixels. The inference step is also decomposed in two distinct parts which are interdependent. The first part consists in computing the approximate posterior distribution, while the second part involves learning the deterministic parameters of the model in θ by maximizing the likelihood of the pixels.

Note that this process leads to an unsupervised learning problem, that is, the model parameters in θ and posterior distribution are both computed based only on the observed HI pixels (i.e., there is no separate training and testing data). In the following, we provide a high-level description of the approach, which is illustrated in Figures 1 and2. Modeling and inference steps are then detailed in Sections IV and V.

Modeling step: The first step is to characterize the dynamical evolution of the EMs and abundances. Under a Markovity assumption, it can be expressed using the following sequence of conditional probability distributions [START_REF] Särkkä | Bayesian filtering and smoothing[END_REF]:

pa t , M t q " p θ pa t , M t |a t´1 , M t´1 q , (2)
where a t " ra J 1,t , . . . , a J N,t s J denotes the abundance maps in lexicographic ordering and M t " M 1,1 , . . . , M N,t ( the collection of EM matrices for all pixels, and pa 0 , M 0 q " p θ pa 0 , M 0 q. We also assume that the abundances and endmembers at time t are statistically independent when conditioned on the their values at time t ´1, that is,

p θ pa t , M t |a t´1 , M t´1 q " p θ pa t |a t´1 qp θ pM t |M t´1 q . (3)
This allows us to model the evolution of a t and M t separately.

The second part of the model represents how the HI pixels are generated from a t and M t , which is given by

y t " p θ py t |a t , M t q , (4)
where y t " ry J 1,t , . . . , y J N,t s J denotes the HI in lexicographic ordering. Note that the pixels y t are assumed to be conditionally independent given a t and M t . These PDFs are defined explicitly in Section IV.

Inference step: This step, which constitutes the solution to the MTHU problem, consists in computing the posterior PDF of the EMs and of the abundances given the observed pixels, which is given by:

p θ pa 1 , M 1 , . . . , a T , M T |y 1 , . . . , y T q (5)
" p θ pa 0 , M 0 q ś T t"1 p θ py t |a t , M t qp θ pa t |a t´1 qp θ pM t |M t´1 q p θ py 1 , . . . , y T q , where the r.h.s. of ( 5) was obtained using the Bayes rule and the factorization in ( 2)-( 4). The PDF in (5) generally does not have a closed form solution [START_REF] Särkkä | Bayesian filtering and smoothing[END_REF]. One efficient solution is to use variational inference based on SGD [START_REF] Kingma | Auto-encoding variational bayes[END_REF], which attempts to find an approximate posterior q P Q within a family of distributions Q that is as close as possible to the true posterior in ( 5). This approximation is often obtained by maximizing a lower bound on the marginal log-likelihood, the so-called evidence lower bound (ELBO):

ELBOpq, θq ď log p θ py 1 , . . . , y T q ,

see Section V for a detailed explanation. Under specific conditions over the model and posterior family Q, such as assuming conditionally Gaussian distributions, the maximization max qPQ ELBOpq, θq can be solved locally using SGD techniques, which are computationally efficient when compared to solutions based on Monte Carlo sampling [START_REF] Kingma | Auto-encoding variational bayes[END_REF].

The flexibility provided by the family of posterior distributions Q is paramount for the performance of the strategy. Recent works have considered neural networks, parameterized by φ, to represent the approximate posterior. Thus, instead of searching for q in a (continuous) family of distributions Q, we search for the parameters φ, such that the parameterized posterior, denoted by q φ , maximizes the ELBO. Thus, the optimization becomes max φ ELBOpq φ , θq. For problems with a temporal Markov structure, RNNs provide a parameterization that, although flexible, is computationally efficient. Moreover, they explicitly explore the temporal structure in the data, having shown excellent performance in various sequence modeling tasks [START_REF] Lipton | A critical review of recurrent neural networks for sequence learning[END_REF]. This will motivate us to use an RNN in the parametrization of our posterior in Section V-C.

Finally, the parameters of the generative model, θ, are also learned within the same framework. The underlying idea is to perform type-II maximum likelihood (ML) estimation [START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF], that is, max θ log p θ py 1 , . . . , y T q. However, since computing p θ py 1 , . . . , y T q is intractable, θ is also computed by maximizing the ELBO w.r.t. θ using SGD. Thus, as ELBOpq φ , θq is maximized w.r.t. q φ , the lower bound in (6) becomes tighter and its maximization w.r.t. θ better approximates ML estimation. Thus, the complete inference problem is formulated as the maximization of the ELBO w.r.t. both the posterior and the model parameters, that is, max θ, φ ELBOpq φ , θq.

IV. PROPOSED MODEL

The modeling step will be divided as follows. First, we develop a mixing model and represent EM variability (which defines p θ py t |a t , M t q). Next, we consider the dynamical behavior of the EMs, and finally of the abundances (which define p θ pM t |M t´1 q and p θ pa t |a t´1 q, respectively).

A. Mixture model with EM variability

As discussed in Section II, devising EM models that combine flexibility to represent complex spectral variability with simplicity of having a small number of parameters is challenging. Flexible models such as the PLMM [START_REF] Thouvenin | Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[END_REF], GLMM [START_REF] Imbiriba | Generalized linear mixing model accounting for endmember variability[END_REF] and ALMM [START_REF] Hong | An augmented linear mixing model to address spectral variability for hyperspectral unmixing[END_REF] have many degrees of freedom and require additional regularization strategies to guarantee physically meaningful solutions, while simpler models such as the ELMM [START_REF] Drumetz | Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability[END_REF] are too restrictive to represent complex spectral variability. An important information which can be used in the design of mixing models accounting for endmember variability is the spectral correlation of EM signatures [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF]. This points to a natural representation of EMs as smooth functions. Although the smoothness of the EMs can be introduced through regularization (see, e.g., [START_REF] Borsoi | Super-resolution for hyperspectral and multispectral image fusion accounting for seasonal spectral variability[END_REF], [START_REF] Halimi | Hyperspectral unmixing in presence of endmember variability, nonlinearity, or mismodeling effects[END_REF]), this leads to highdimensional and potentially costly HU solutions. On the other hand, a more efficient and interpretable model can be obtained by directly parametrizing smoothness using properly selected basis functions [START_REF] Halimi | Fast hyperspectral unmixing in presence of nonlinearity or mismodelling effects[END_REF].

In this work, we consider an EM model inspired by the GLMM [START_REF] Imbiriba | Generalized linear mixing model accounting for endmember variability[END_REF], which represents spectral variability using a multiplicative scaling of reference EM spectra that vary for each band, endmember and pixel. However, instead of using regularizations we propose to constrain the scaling factors to be linear combinations of spectrally smooth functions. The resulting model, which we call Smooth GLMM (SGLMM), represents each observed HI pixel y n,t as follows:

y n,t " `M 0 d p1 `D r Ψ n,t q looooooooooooomooooooooooooon M n,t a n,t `rn,t , (7)
where d represents the Hadamard (elementwise) product, 1 is an L ˆP matrix of ones, M 0 P R LˆP a set of reference or average EM signatures, matrix D P R LˆK contains K spectrally smooth basis vectors as its columns, and r Ψ n,t P R KˆP contains the low-dimensional coefficients that parameterize the variability of each EM. Vector r n,t P R L denotes zero-mean additive Gaussian noise.

It is instructive to analyze how the variability of the endmembers M n,t is introduced in the model [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF]. The EM matrix M n,t is formed by scaling the reference EMs in M 0 bandwise by matrix 1 `D r Ψ n,t P R LˆP . This model is similar to the GLMM, the difference being in the structure of this multiplicative scaling matrix. First, note that when r Ψ n,t « 0, the term D r Ψ n,t is also small and the scaling factors will be close to 1, meaning that M n,t « M 0 , i.e., the spectral variability is small. Thus, the amount of EM variability depends directly on the amplitude of the elements of r Ψ n,t . Second, matrix D r Ψ n,t represents a perturbation over the constant scaling 1, and its properties depend directly on the choice of D. Thus, by properly selecting D we can constrain D r Ψ n,t to represent smooth functions with few parameters, leading to smooth spectral variations in M n,t . Following an idea used in [START_REF] Halimi | Fast hyperspectral unmixing in presence of nonlinearity or mismodelling effects[END_REF] for robust HU with smooth additive residual terms, we select the columns of D as the first K rows of the discrete cosine transform (DCT) matrix.

The SGLMM models endmember variability using KP parameters. The number of basis functions K gives a tradeoff between existing models in the literature: when K " 1, D will contain only a constant vector and the model becomes equivalent to the ELMM [START_REF] Drumetz | Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability[END_REF], whereas for K " L it has the same flexibility as the PLMM [START_REF] Thouvenin | Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[END_REF] and GLMM [START_REF] Imbiriba | Generalized linear mixing model accounting for endmember variability[END_REF]. Values of K ! L should give the SGLMM sufficient flexibility to represent smooth spectral variability accurately.

Representing the spectral variability parameters in vectorized form as ψ n,t " vecp r Ψ n,t q and assuming the noise r n,t to be independent for each pixel, the PDFs in the generative model ( 4) can be rewritten equivalently in terms of ψ n,t as

p θ py t |a t , ψ t q " N ź n"1 p θ py n,t |a n,t , ψ n,t q , (8)
where

p θ py n,t |a n,t , ψ n,t q " N ´`M 0 d p1 `D vec ´1pψ n,t qq ˘an,t , σ 2 r I ¯, (9)
in which σ r P R ˚is the standard deviation of the measurement noise, which is assumed to be independent and identically distributed for different bands.

B. Dynamical model for the EM scaling parameters

In this work, we consider M 0 to be a deterministic parameter of the model and estimated from the observed HIs using an approximate ML framework (i.e., M 0 P θ). This means that the EM matrix for each t and n, M n,t " M 0 d p1 `D r Ψ n,t q, is a deterministic function of the lowerdimensional vector of scaling factors ψ n,t . Thus, we can substitute the problem of estimating the very high-dimensional p θ pM t |M t´1 q by the problem of estimating p θ pψ t |ψ t´1 q, with ψ t " " ψ J 1,t , . . . , ψ J N,t ‰ J . Since vector ψ t is still high dimensional, we consider an independence assumption on the time evolution between different pixels:

p θ pψ t |ψ t´1 q " N ź n"1 p θ pψ n,t |ψ n,t´1 q , ( 10
)
for t ě 1, where the prior PDF for time instant t " 0 will be specified later in Section IV-D. We consider a Gaussian distribution to represent the evolution of ψ n,t :

p θ pψ n,t |ψ n,t´1 q " N `ψn,t´1 , σ 2 ψ I P K ˘, (11)
where σ ψ P R ˚is the distribution standard deviation, which controls its uncertainty and is assumed to be isotropic. Note that the evolution of ψ n,t is not assumed to be affected by abrupt changes, which leads us to consider σ ψ constant. This assumption is motivated from the fact that the reflectance of materials are primarily influenced by their physico-chemical composition (e.g., particle size and roughness in packed particle spectra [START_REF] Hapke | Bidirectional reflectance spectroscopy, 1, Theory[END_REF], or biophysical parameters in leaf spectra [START_REF] Jacquemoud | Leaf optical properties: A state of the art[END_REF]), which we assume to change smoothly at fine time scales.

C. Abundances model

In order to represent the abundances dynamical behavior, we first assume their time evolution to be independent for different pixels in order to make the problem tractable, that is,

p θ pa t |a t´1 q " N ź n"1 p θ pa n,t |a n,t´1 q , ( 12
)
for t ě 1, where the prior PDF for time instant t " 0 will be specified later in Section IV-D. To represent the time evolution at each pixel, we consider a Dirichlet distribution.

The Dirichlet is a natural choice of distribution to model the abundances as it enforces the physical constraints that the elements of a t should be nonnegativity and sum to one [START_REF] Halimi | Unsupervised unmixing of hyperspectral images accounting for endmember variability[END_REF], [START_REF] Eches | Adaptive Markov random fields for joint unmixing and segmentation of hyperspectral images[END_REF]. The transition PDF is then given by p θ pa n,t |a n,t´1 q " Dirpα n,t q , (13) where α n,t P R P `denotes the concentration parameters, which are a function of the abundances at the previous time instant, a n,t´1 (i.e., the parameters of p θ pa n,t |a n,t´1 q, α n,t , are a function of the conditioning variable). Note that the uncertainty of the abundances predictions is represented implicitly in α n,t , where small concentration values yield low uncertainty and temporally smooth transition, whereas large values lead to higher uncertainty, allowing for more changes.

However, the Dirichlet distribution can make inference difficult. One workaround consists of using, e.g., Laplace's method, which approximates the Dirichlet distribution Dirpα n,t q by a Gaussian with mean and inverse covariance equal to the mode of the original distribution and the Hessian of its negative logarithm, respectively [START_REF] Barber | Bayesian Reasoning and Machine Learning[END_REF]. However, since the Dirichlet distribution is supported at the simplex, this approximation can be inaccurate. To overcome this problem, MacKay [START_REF] Mackay | Choice of basis for Laplace approximation[END_REF] proposed to perform this approximation in the so-called softmax basis, which consists in a mapping π ´1 : a n,t Þ Ñ c n,t from the unity simplex to R P , where π is the softmax function:

π ´1pa n,t q " c n,t , (14)
π i pc n,t q " exppc n,t,i q ř j exppc n,t,j q , i P t1, . . . , P u ,

with π i , c n,t,i being the i-th positions of π and c n,t . This approximation is very accurate and has been used in several works to facilitate statistical inference [START_REF] Hennig | Kernel topic models[END_REF], [START_REF] Srivastava | Autoencoding variational inference for topic models[END_REF]. Thus, replacing a n,t by the softmax parameters c n,t , we achieve the following alternative representation of (13) [START_REF] Mackay | Choice of basis for Laplace approximation[END_REF]:

p θ pc n,t |c n,t´1 q " Γp ř P i"1 α n,t,i q ś P i"1 Γpα n,t,i q P ź i"1 π i pc n,t q αn,t,i gp1 J c n,t q , ( 16
)
where α n,t,i is the i-th position of α n,t , Γ denotes the Gamma function, and g is an arbitrary distribution used to constrain an extra degree of freedom (since the Dirichlet has only P ´1 degrees of freedom), selected as gpxq9 expp´ 2 x 2 q for mathematical convenience [START_REF] Mackay | Choice of basis for Laplace approximation[END_REF]. The Gaussian approximation of this distribution is then given by p θ pc n,t |c n,t´1 q « N pµ n,t , Σ n,t q [26], [START_REF] Hennig | Kernel topic models[END_REF], with µ n,t given by µ n,t,i " logpα n,t,i q ´1

P P ÿ "1 logpα n,t, q , (17)
where µ n,t,i is the i-th position of µ n,t and Σ n,t is the negative Hessian of ( 16) at c n,t " µ n,t . Note that the mean and covariance µ n,t and Σ n,t are a function of α n,t and, consequently, depend implicitly on c n,t´1 . Therefore, we can approximate the transition PDF ( 13) by a Gaussian one on the softmax basis. Note that the relationship between the parameters of both models, that is, between α n,t and µ n,t and Σ n,t , is nonlinear and burdensome to compute. However, by working on the softmax basis we do not need to specify p θ pa n,t |a n,t´1 q explicitly in [START_REF] Goenaga | Unmixing analysis of a time series of hyperion images over the guánica dry forest in puerto rico[END_REF]. Instead, we can directly define the Gaussian transition for p θ pc n,t |c n,t´1 q, which is mathematically more convenient. This implicitly defines a transition probability p θ pa n,t |a n,t´1 q by mapping c n,t " p θ pc n,t |c n,t´1 q into the simplex, approximating a Dirichlet distribution. Thus, we consider the following model:

p θ pc n,t |c n,t´1 q " N `cn,t´1 , σ 2 a pc n,t´1 qI P ˘, (18)
where I P is a P ˆP identity matrix. Note that, for simplicity, the covariance matrix in [START_REF] Palsson | Blind hyperspectral unmixing using autoencoders: A critical comparison[END_REF] was constrained to be isotropic, and is scaled by σ 2 a pc n,t´1 q. The function σ a : R P Ñ R computes the standard deviation of each element of c n,t " p θ pc n,t |c n,t´1 q as a function of c n,t´1 . Thus, it directly influences the amount of change in the abundances: the larger σ a pc n,t´1 q, the larger the changes we expect to observe between c n,t´1 and c n,t . This function, which is part of the generative model, will be learned during inference using a maximum likelihood approach. It will be parameterized using a fully connected neural network with R σa layers, where each hidden layer has P neurons and uses the ReLU activation function, and the output layer maps to a scalar and uses an exponential activation function to ensure the output is positive.

D. The complete model

To finish the model derivation, we need to define the initial PDFs at time t " 0, which under the new parametrization of the abundances and endmembers, which we assume to be pixelwise independent Gaussian distributions, given by:

p θ pc 0 , ψ 0 q " N ź n"1 p θ pc n,0 qp θ pψ n,0 q , ( 19
)
p θ pc n,0 q " N `νc 0 , diagpγ c 0 q 2 ˘, (20)
p θ pψ n,0 q " N `νψ 0 , diagpγ ψ 0 q 2 ˘, (21)
for all n " 1, . . . , N , where the means ν c 0 , ν ψ 0 and diagonal covariance parameters, γ c 0 , γ ψ 0 are constant and shared among all pixels in order to reduce the amount of parameters in the model. Finally, the measurement model ( 9) can be written using the softmax abundance reparametrization as:

p θ py n,t |c n,t , ψ n,t q " N ´`M 0 d p1 `D vec ´1pψ n,t qq ˘πpc n,t q, σ 2 r I ¯, ( 22
)
where π is the softmax function.

The final dynamical model is then given by equations [START_REF] Wang | Nonlinear unmixing of hyperspectral data via deep autoencoder networks[END_REF], (21) (initial PDFs), [START_REF] Palsson | Blind hyperspectral unmixing using autoencoders: A critical comparison[END_REF], [START_REF]Multi-temporal hyperspectral mixture analysis and feature selection invasive species mapping in rainforests[END_REF] (the dynamical model) and ( 22) (the measurement model). Finally, we denote the parameters of the model which will be estimated from the data using approximate ML inference by θ " tM 0 , σ r , σ a , ν c 0 , ν ψ 0 , γ c 0 , γ ψ 0 u. An illustrative diagram of the proposed generative model can be seen in Figure 3.

V. VARIATIONAL INFERENCE WITH RNNS FOR HU

In this section, we will present the proposed solution to the inference step, referred to as ReDSUNN. Considering the parametrization of the abundances and of endmember variability derived in the previous section, this task, which consists in performing MTHU, becomes that of approximating the posterior distribution:

p θ `c1 , ψ 1 , . . . , c T , ψ T ˇˇy 1 , . . . , y T ˘, (23)
where c t " " c J 1,t , . . . , c J N,t ‰ J . First, let us denote with an underline the collection of variables at all time instants: y " ty 1 , . . . , y T u , (24) y n " ty n,1 , . . . , y n,T u ,

c " tc 0 , . . . , c T u , (26) ψ " tψ 0 , . . . , ψ T u .

(

) 27
As discussed in Section III, due to the nonlinearity in the model caused by the interaction between the abundances and the variability scaling factors, and the potentially high dimensionality of these variables, it is not possible to compute the posterior distribution [START_REF] Barber | Bayesian Reasoning and Machine Learning[END_REF] in closed form. In this work, we adopt a deep variational inference framework: we consider a parametric surrogate distribution q φ pc, ψ|yq from a sufficiently flexible family with parameters φ, and learn its parameters by minimizing the Kullback-Leibler (KL) divergence between q φ pc, ψ|yq and the true posterior p θ pc, ψ|yq:

KL `pθ pc, ψ|yq › › q φ pc, ψ|yq ˘" log p θ pyq `Eq φ pc,ψ|yq tlog q φ pc, ψ|yqu ´Eq φ pc,ψ|yq tlog p θ pc, ψ, yqu .

Since the KL divergence is nonnegative and log p θ pyq is a constant, the above expression can be equivalently minimized by maximizing a lower bound to the data likelihood formed by the last two terms in the right hand side of the expression, which is the so-called ELBO [START_REF] Barber | Bayesian Reasoning and Machine Learning[END_REF], [START_REF] Kingma | Auto-encoding variational bayes[END_REF]:

log p θ pyq ě E q φ pc,ψ|yq tlog p θ pc, ψ, yqu ´Eq φ pc,ψ|yq tlog q φ pc, ψ|yqu .

Recent advances in variational deep learning such as in variational autoencoders has made it possible to devise efficient algorithms to maximize (29) when the PDFs are possibly parameterized by deep neural network using, e.g., stochastic backpropagation algorithms [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]. Furthermore, the conditional independence assumptions of the model in the previous section can be exploited to further simplify the inference problem. In the following, we factorize the ELBO both in the temporal as well as in the pixel dimensions. Note that, as discussed in Section III, (29) will be optimized both with respect to the parameters of the posterior distribution φ, but also with respect to the parameters of the generative model p θ pc, ψ, yq in the set θ in order to estimate them by approximate ML inference. In the following, we will use the Markov and pixelwise conditional independence assumptions of the generative model in Section IV to factorize q φ pc, ψ|yq and simplify the solution to the inference problem.

A. Factorizing the posterior distribution

Factorizing the posterior distribution in time: Various kinds of parametrizations of the distribution q φ pc, ψ|yq have been proposed. One of the simplest is to consider a mean field assumption [START_REF] Krishnan | Structured inference networks for nonlinear state space models[END_REF], which assumes that tc t , ψ t u and tc t 1 , ψ t 1 u are conditionally independent given y, for all t ‰ t 1 . However, this disregards the temporal structure of the data. A more suitable factorization can be obtained by noting that the Markov property of the model can be used to show that the true posterior factorizes as

p θ pc, ψ|yq " p θ pc 0 , ψ 0 |yq T ź t"1 p θ pc t , ψ t |c t´1 , ψ t´1 , yq .
Incorporating this assumption into the variational approximation q φ pa, ψ|yq leads to a similar factorization:

q φ pc, ψ|yq " q φ pc 0 , ψ 0 |yq T ź t"1 q φ pc t , ψ t |c t´1 , ψ t´1 , yq , (30)
which preserves the temporal dependency of the model.

Factorizing the posterior distribution in pixels: The vectors c t , ψ t in (30) contain the abundances and variability coefficients for all image pixels, and are thus of very high dimension. Therefore, additional simplifications are necessary in order to make inference tractable. One important property is that in the model derived in Section IV, the initial, transition, and measurement PDFs (equations ( 20), ( 21), ( 18), [START_REF]Multi-temporal hyperspectral mixture analysis and feature selection invasive species mapping in rainforests[END_REF] and ( 22)) can be factorized among the different image pixels. Thus, the inference process can be factorized at the pixel level, which leads to the following form for the posterior distribution:

q φ pc t , ψ t |c t´1 , ψ t´1 , yq " N ź n"1 q φ pc n,t , ψ n,t |c n,t´1 , ψ n,t´1 , y n q , ( 31
)
for t ě 1, and similarly for the initial PDF at t " 0:

q φ pc 0 , ψ 0 |yq " N ź n"1 q φ pc n,0 , ψ n,0 |y n q . ( 32
)
Although this factorization does not directly consider spatial correlation between different pixels, which has been found to be a useful source of prior information in HU [START_REF] Eches | Enhancing hyperspectral image unmixing with spatial correlations[END_REF], [START_REF] Halimi | Unsupervised unmixing of hyperspectral images accounting for endmember variability[END_REF], it allows us to work with pixelwise variational posterior PDFs (i.e., the r.h.s. of ( 31) and ( 32)) which have a much lower dimension, thus reducing the computational burden associated with the inference step. To some extent, spatial information can still be introduced indirectly by constraining the parametrization of the variational posterior distributions among different pixels, which will be explained in the rest of this subsection. The incorporation of spatial information directly through the probabilistic model will be investigated in a future work.

Parameterizing the posterior: A key aspect of the model is how to parameterize the posterior PDFs of the different pixels in the r.h.s. of ( 31) and [START_REF] Roberts | Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[END_REF]. First, variational inference implies selecting a parametric family of distributions from which to select q φ , which directly impact the results. Note that the true posterior in (5) might have a complex form and be possibly multimodal, however, its form is not known in advance. Thus, as in recent works in deep variational inference (see, e.g., [START_REF] Kingma | Auto-encoding variational bayes[END_REF]) we considered a Gaussian family for q φ since this will simplify the maximization of the ELBO considerably (through, e.g., the reparametrization trick and closed form expressions for KL divergences), leading to important computation savings. Thus, it can be expressed as:

q φ `cn,t , ψ n,t ˇˇc n,t´1 , ψ n,t´1 , y n " N `µc,ψ φ pΥ n,t q, diagpσ c,ψ φ pΥ n,t qq 2 ˘, (33)
where Υ n,t " c n,t´1 , ψ n,t´1 , y n ( and µ c,ψ φ and σ c,ψ φ are functions (e.g., neural networks parameterized by φ) which compute the parameters of the posterior distribution, mapping the data tc n,t´1 , ψ n,t´1 , y n u to the mean and the square root of the diagonal covariance matrix of the Gaussian posterior, respectively. For convenience of notation, we decompose µ c,ψ φ and σ c,ψ φ into two functions:

µ c,ψ φ " " µ c φ µ ψ φ  , σ c,ψ φ " " σ c φ σ ψ φ  . (34)
Note that functions µ c φ and σ c φ compute the mean and square root of the diagonal covariance matrix of the q φ `cn,t ˇˇc n,t´1 , ψ n,t´1 , y n ˘, while µ ψ φ and σ ψ φ compute the mean and square root of the diagonal covariance matrix of q φ `ψn,t ˇˇc n,t´1 , ψ n,t´1 , y n ˘. A Gaussian parametrization is also used for the posterior distribution of the initial PDF:

q φ pc 0 , ψ 0 |yq " N `ζc,ψ , diagpξ c,ψ q 2 ˘, (35)
where a fixed distribution was used for all pixels, with ζ c,ψ and ξ c,ψ being the mean and the diagonal of the square root of the covariance matrix, respectively. An important observation is that we consider a shared parametrization, where the posterior in [START_REF] Dudley | A multi-temporal spectral library approach for mapping vegetation species across spatial and temporal phenological gradients[END_REF] and [START_REF]Real-time spatiotemporal spectral unmixing of MODIS images[END_REF] has the same form for all pixels. More precisely, this means that the same functions µ c,ψ φ and σ c,ψ φ are used to compute the posterior mean and covariance for every HI pixel, given the input data Υ n,t . This is an important characteristic of the method, since it allows information from different pixels (i.e., from the whole image) to be leveraged jointly in the estimation of the model and, consequently, of the abundances and variability coefficients in each pixel, c n,t , ψ n,t , n " 1, . . . , N .

B. Factorizing the ELBO cost function

Using the simplifications derived in the previous subsection, in the following we will rewrite the ELBO cost function [START_REF] Sigurdsson | Sparse distributed multitemporal hyperspectral unmixing[END_REF] in terms of the factorized model.

Factorizing the ELBO temporally: Using the factorization [START_REF] Thouvenin | A hierarchical bayesian model accounting for endmember variability and abrupt spectral changes to unmix multitemporal hyperspectral images[END_REF] and the Markovity of the model, the lower bound in ( 29) can be written as [START_REF] Krishnan | Structured inference networks for nonlinear state space models[END_REF]:

log p θ pyq ě Lpθ, φ, yq " T ÿ t"1 E q φ pct,ψ t |yq log p θ py t |c t , ψ t q ( ´KL `qφ pc 0 , ψ 0 |yq › › p θ pa 0 , ψ 0 q ˘´T ÿ t"1 E q φ pct´1,ψ t´1 |yq ! KL `qφ pc t , ψ t |c t´1 , ψ t´1 , yq › › p θ pc t , ψ t |c t´1 , ψ t´1 q ˘) . ( 36
)
Factorizing at pixel level: Using the pixelwise factorization of the generative and posterior PDFs discussed in Section IV, we can simplify each term of [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF]. To this end, we use the fact that KLpppx 1 , x 2 q}qpx 1 , x 2 qq " KLpppx 1 q}qpx 1 qq KLpppx 2 q}qpx 2 qq when both ppx 1 , x 2 q " ppx 1 qppx 2 q and qpx 1 , x 2 q " qpx 1 qqpx 2 q are independent, and the fact that E ppx1,x2q tf px 1 qu " E ppx1q tf px 1 qu. We proceed to analyse each term of (36) in the following.

First term: Using the factorization of the measurement model in [START_REF] Borsoi | Kalman filtering and expectation maximization for multitemporal spectral unmixing[END_REF], the first term in the r.h.s. of [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] becomes:

E q φ pct,ψ t |yq log p θ py t |c t , ψ t q ( " N ÿ n"1 E q φ pcn,t,ψ n,t |y n q log p θ py n,t |c n,t , ψ n,t q ( . (37)
Second term: Using the pixelwise independence of the initial PDF in the generative model [START_REF] Zhou | ADMM-Based hyperspectral unmixing networks for abundance and endmember estimation[END_REF] and in the posterior [START_REF] Roberts | Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[END_REF], the KL divergence can be written as:

KL `qφ pc 0 , ψ 0 |yq › › p θ pc 0 , ψ 0 q " N ÿ n"1 KL `qφ pc n,0 , ψ n,0 |y n q › › p θ pc n,0 , ψ n,0 q ˘. ( 38
)
Third term: Using the pixelwise independence of the posterior [START_REF] Zhuo | A novel spectral-temporal Bayesian unmixing algorithm with spatial prior for Sentinel-2 time series[END_REF] and of the predictive PDFs ( 12), [START_REF] Somers | Invasive species mapping in hawaiian rainforests using multi-temporal hyperion spaceborne imaging spectroscopy[END_REF], this term can be written as:

E q φ pct´1,ψ t´1 |yq KL qφ pc t , ψ t |c t´1 , ψ t´1 , yq › › p θ pc t , ψ t |c t´1 , ψ t´1 q ˘( " N ÿ n"1 E q φ pcn,t´1,ψ n,t´1 |y n q KL `(39) q φ pc n,t , ψ n,t |c n,t´1 , ψ n,t´1 , y n q › › p θ pc n,t , ψ n,t |c t´1 , ψ t´1 q ˘(.
Combining these results, we can write the cost function Lpθ, φ, yq as in equation ( 40) (depicted on top of the next page). Details on the computation of the log-likelihood and KL divergences can be found in Appendix A.

C. An RNN-based implementation

A key question is how to define the functions µ c φ , µ ψ φ , σ c φ and σ ψ φ in ( 33) and [START_REF] Wang | Spatio-temporal spectral unmixing of time-series images[END_REF], which parameterize the approximate posterior distribution. On the one hand, these have to be flexible to be able to approximate the true posterior, which cannot be written in the form of a simple and well-known Bidirectional RNN Fig. 4: Diagram of the proposed network implementing the posteriors q φ pcn,t, ψ n,t |c n,t´1 , ψ n,t´1 , y n q and q φ pc 0 , ψ 0 |yq.

distribution. On the other hand, it is important to incorporate information from a physical modeling of the problem to make inference process more efficient, interpretable and stable. Thus, we will parameterize the variational posterior distribution using a lightweight RNN and, whenever possible, leveraging physically motivated models. First, a bidirectional RNN is used to compute a set of feature-based representations, denoted by h n,t P R H , from the image pixel sequence y n . In particular, we compute h n,t by combining the hidden states learned by two LSTMs [START_REF] Hochreiter | Long short-term memory[END_REF]:

h forw n,t " LSTM forw φ ph forw n,t´1 , y n,t q, t " 1, . . . , T , (41)
h back n,t " LSTM back φ ph back n,t`1 , y n,t q, t " T, . . . , 1 , (42)
h n,t " 1 2 `hforw n,t `hback n,t ˘, (43)
for n " 1, . . . , N , where LSTM forw φ and LSTM back φ denote two LSTMs which process the data forward and backwards in time, respectively; their hidden state representation being given by h forw n,t and h back n,t . We choose LSTMs due to their excellent performance in various sequence modeling tasks [START_REF] Lipton | A critical review of recurrent neural networks for sequence learning[END_REF]. Moreover, a bidirectional RNN (i.e., two LSTMs) is used because at every time instant the posterior in [START_REF] Dudley | A multi-temporal spectral library approach for mapping vegetation species across spatial and temporal phenological gradients[END_REF] depends on the HI pixels at all time instants, y n , whereas the LSTMs in [START_REF] Halimi | Unsupervised unmixing of hyperspectral images accounting for endmember variability[END_REF] and [START_REF] Zhou | A gaussian mixture model representation of endmember variability in hyperspectral unmixing[END_REF] depend only on past and future data, respectively.

The dimension of the RNN representation is selected as H " pK `1qP (i.e., the dimension of the state vector). The gating units of the LSTMs use the sigmoid nonlinearity, while the input and hidden state units use the uses the hyperbolic tangent nonlinearity. Note that the parameters of these LSTMs will also be learned during inference by SDG using backpropagation through time [START_REF] Lipton | A critical review of recurrent neural networks for sequence learning[END_REF].

We now use the representation h n,t to parameterize µ c φ , µ ψ φ , σ c φ and σ ψ φ . To introduce physical knowledge, we follow the general idea of using hybrid models [START_REF] Imbiriba | Hybrid Neural Network Augmented Physics-based Models for Nonlinear Filtering[END_REF], [START_REF] Li | Model-based deep autoencoder networks for nonlinear hyperspectral unmixing[END_REF], in which an approximate model is complemented by a learnable component (in this case derived from the RNN). In particular, for the posterior mean of the abundances, µ c φ , we construct an approximate model by assuming that 1) a least squares solution provides a crude abundance estimate, 2) the abundances are temporally smooth but may undergo sudden changes, and 3) abrupt abundance changes lead to abrupt changes in the pixels. For the variability parameters, µ ψ φ , we consider it to Lpθ, φ,yq "

T ÿ t"1 N ÿ n"1 E q φ pcn,t,ψ n,t |y n q log p θ py n,t |c n,t , ψ n,t q ( ´N ÿ n"1 KL `qφ pc n,0 , ψ n,0 |y n q › › p θ pc n,0 , ψ n,0 q T ÿ t"1 N ÿ n"1 E q φ pcn,t´1,ψ n,t´1 |y n q ! KL `qφ pc n,t , ψ n,t |c n,t´1 , ψ n,t´1 , y n q › › p θ pc n,t , ψ n,t |c n,t´1 , ψ n,t´1 q ˘) . (40)
be temporally smooth. For the standard deviations, σ c φ and σ ψ φ , we don't have a good physical model; thus, we use a purely non-parametric representation. In the following, we define each of these functions explicitly; an illustrative diagram can be seen in Figure 4.

Considering the RNN features h n,t , the abundance posterior means µ c φ is then parameterized as:

µ c φ pΥ n,t q " π ´1´α 1 p1 ´un,t qπpc n,t´1 q `α2 u n,t `M : n,t´1 y n,t `Wc h n,t ˘¯, (44)
where M n,t´1 " M 0 d p1 `D vec ´1pψ n,t´1 qq is the predicted EM matrix at pixel n and time t ´1, α 1 and α 2 are trainable, real-valued weighting coefficients, and W c P R P ˆpK`1qP is a trainable matrix that maps the hidden RNN representations to the abundances in the softmax basis. The scalar coefficient u n,t P r0, 1s, defined as u n,t " 1 2P }spM : n,t´1 y n,t q ´πpc n,t´1 q} 1 , measures the difference between the predicted abundances πpc n,t´1 q and a crude estimation of the current abundances at time t, given by spM :

n,t´1 y n,t q, where the fixed function sp¨q projects the linear regression solution M : n,t´1 y n,t to the unit simplex. Thus, u n,t works as a crude abrupt change detector.

The parametrization (44) can be seen as a weighted combination of three terms: the abundances at the previous time instant, a crude abundance estimate at time t computed by linear regression, and a non-parametric term depending on h n,t . The balance between them depends on the trainable weights and on the change detector u n,t . When there are no changes, u n,t is small, which gives a higher contribution to the predicted abundances c n,t´1 in [START_REF] Thouvenin | Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[END_REF]. On the other hand, if there is an abrupt change, u n,t is large, giving a higher contribution to the sum of the last two terms in [START_REF] Thouvenin | Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[END_REF], which is a linear regression-based abundance estimate augmented by a non-parametric RNN-based representation. This parametrization is particularly relevant since the generative model [START_REF] Palsson | Blind hyperspectral unmixing using autoencoders: A critical comparison[END_REF] does not explicitly represent abrupt changes.

For the function µ ψ φ , to leverage temporal smoothness we consider a weighted linear combination of the variability coefficients at the previous time instant ψ n,t´1 and a linear mapping of the hidden RNN representation:

µ ψ φ pΥ n,t q " βψ n,t´1 `Wψ h n,t , (45)
where β is a real-valued weight, and W ψ P R KP ˆpK`1qP is a matrix that computes the variability coefficients' innovation from the RNN representation h n,t , both of which are trainable. Note that by not considering abrupt changes to occur in ψ t we obtain a simpler model compared to µ c φ .

The standard deviations σ c φ and σ ψ φ are computed based on a fully non-parametric model, which is given as linear mappings of the RNN representations h n,t :

σ c φ pΥ n,t q " exp `V c h n,t ˘, (46)
σ ψ φ pΥ n,t q " exp `V ψ h n,t ˘, (47)
where V c P R P ˆpK`1qP and V ψ P R KP ˆpK`1qP are the transformation matrices, and the exponential function is applied elementwise in order to ensure nonnegativity of the standard deviations. The parameters of the approximate posterior are finally denoted by φ "

ζ c,ψ , ξ c,ψ , LSTM forw φ , LSTM back φ , α 1 , α 2 , β, W c , V c , W ψ , V ψ ( .
Note that all parameters in φ will be learned using SGD.

Approximating the expectations and optimization: To optimize (40) using stochastic backpropagation, it is necessary to estimate gradients of expectations whose distribution depend on θ and φ, which are the parameters to be optimized. To address this issue, we consider the reparametrization trick, which provides low-variance gradient estimates [START_REF] Kingma | Auto-encoding variational bayes[END_REF]. This is performed by writing the random variables inside the expectations as deterministic functions of a random variable that does not depend on φ. In general, for a distribution q φ pxq and function f , this can be formulated as E q φ pxq tf pxqu " E pp q tf pgp qqu, where g is a function such that x and gp q have the same distribution, and pp q does not depend on φ. Applying this to the expectations in [START_REF] Uezato | Hyperspectral unmixing with spectral variability using adaptive bundles and double sparsity[END_REF] and considering that the posterior in our model is Gaussian, this is achieved as:

" a J n,t , ψ J n,t ‰ J " µ c,ψ φ pΥ n,t q `σc,ψ φ pΥ n,t q d , (48)
for t " 1, . . . , T , where " N p0, Iq. It can be verified that the random variables in [START_REF] Borsoi | A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF] are sampled according to the distribution q φ pc n,t , ψ n,t |y n q. Thus, by using this reparametrization, the expectations in ( 40) can be rewritten in terms of expectations of , which we subsequently approximate using a one-sample Monte Carlo estimate and denote by p Lpθ, φ, yq.

The approximated cost function p

Lpθ, φ, yq is then optimized with respect to both θ and φ (i.e., the parameters of the generative model and of the variational posterior) using the Adam stochastic optimization method [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]. We used a learning rate of 0.001 and a batch size of 128. Training was performed for 30 epochs. The full MTSU process performed by ReDSUNN is summarized in Algorithm 1.

Since the cost function is non-convex, the initialization of the parameters can have an important impact on the solution. The parameters of the neural networks σ a , LSTM forw φ , LSTM back φ , and the matrices V c and V ψ are initialized randomly using Glorot initialization [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. W c and W ψ are initialized with zeros, and β " α 1 " α 2 " 1. M 0 was initialized using the vertex component analysis (VCA) algorithm [START_REF] Nascimento | Vertex Component Analysis: A fast algorithm to unmix hyperspectral data[END_REF], and σ r " 0.0001 (corresponding to an SNR of about 35dB Algorithm 1: ReDSUNN Input : HIs y 1 , . . . , y T , hyperparameters P , K and σ ψ . 1 Initialize θ 0 and φ 0 as described in Section V-C ; 2 Use Adam [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] to maximize p Lpθ, φ, yq w.r.t. both φ and θ ; 3 for n " 1, . . . , N do 4 Compute p c n,0 , p ψ n,0 as the means of q φ pc n,0 , ψ n,0 |y n q using (35); 5 Compute the h n,1 , . . . , h n,T using ( 41), ( 42) and ( 43 for spectra with standard deviation 0.5). For the parameters of the initial prior and variational posterior PDFs, we initialized the means ν a 0 , ν ψ 0 and ζ c,ψ with zeros, and the variances γ a 0 , γ ψ 0 and ξ c,ψ with ones, making the initial PDFs standard Gaussians.

D. Model complexity and comparisons

We now summarize the parameters of the generative model, of the variational posterior, and their dimensionality (i.e., the number of parameters that have to be inferred). This can be seen in Table I. To compute the number of parameters corresponding to the LSTMs, we note that each LSTM has four input-hidden weight matrices, four hidden-hidden weight matrices, and four biases (where the input is of size L, and the hidden state of size pK `1qP ). It is instructive to compare the amount of parameters to other methods in the literature. By using a Markovity assumption, a shared posterior distribution for all pixels, and an RNN posterior parametrization, the amount of parameters to be learned by ReDSUNN in Table I does not scale with either N or T , differently from previous methods such as OU [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF] or the HBUN [START_REF] Liu | Bayesian unmixing of hyperspectral image sequence with composite priors for abundance and endmember variability[END_REF].

M 0 LP σr 1 σa P pP `1qRσ a ν a 0 , ν ψ 0 , γ a 0 , γ ψ 0 2pK `1qP
Variational posterior (φ)

LSTM forw φ , LSTM back φ 8pK `1qP `pK `1qP `L `1α 1 , α 2 , β 3 Wc, V c 2pP 2 pK `1qq W ψ , V ψ 2pKP 2 pK `1qq ζ c,ψ , ξ c,ψ 2pK `1qP

VI. RESULTS

The performance of the proposed ReDSUNN algorithm is evaluated using simulations with synthetic and real data. We compare our method with the fully constrained least squares (FCLS), online unmixing (OU) [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF], HBUN [START_REF] Liu | Bayesian unmixing of hyperspectral image sequence with composite priors for abundance and endmember variability[END_REF], and with a Kalman filter and expectation maximization-based strategy (referred to simply as Kalman) [START_REF] Borsoi | Kalman filtering and expectation maximization for multitemporal spectral unmixing[END_REF]. The EMs used by FCLS were extracted by the VCA algorithm at each time instant [START_REF] Nascimento | Vertex Component Analysis: A fast algorithm to unmix hyperspectral data[END_REF].

The reference EMs required by the Kalman method, and the initialization for the EMs in OU, HBUN and for the proposed method were all set with the signatures obtained by applying VCA to the matrix ry 1,1 , . . . , y n,t . . . , y N,T s P R LˆN T formed by concatenating the HI pixels for all time instants. The abundances and EM scaling factors estimated by ReDSUNN are set according to Algorithm 1. The hyperparameters of all algorithms were adjusted so as to obtain high abundance reconstruction performance. For ReDSUNN, parameters K and σ ψ (which are not optimized) were searched within the ranges K P t1, . . . , 10u and σ ψ P t10 ´5, . . . , 0.1, 1u, and R σa " 2 layers were used to parameterize function σ a p¨q in [START_REF] Palsson | Blind hyperspectral unmixing using autoencoders: A critical comparison[END_REF]. For the other algorithms, their parameters were selected in the ranges indicated in their original publications. The proposed method was implemented in Pytorch (codes will be available at https://github.com/ricardoborsoi/ReDSUNN). The remaining methods were implemented in Matlab (codes were provided by the original authors). All experiments were run in a desktop computer with an Intel Xeon TM W-2104 CPU with four 3.2GHz cores and 24GB of RAM. No GPU was used in the simulations. ReDSUNN, OU and Kalman used parallelization in their implementations.

The quantitative performance of the algorithms was evaluated using the average normalized mean squared error (NRMSE), between the abundances, EMs, and reconstructed HIs, which are computed as

NRMSE A " `1 T ř T t"1 ř N n"1 }a n,t ´p a n,t } 2 L }a t } 2 ˘1{2 , NRMSE M " `1 N T ř T t"1 ř N n"1 }M n,t ´x M n,t } 2 F L }M n,t } 2 F ˘1{2 , and NRMSE Y " `1 T ř T t"1 ř N n"1 › › y n,t ´x M n,t p a t › › 2 L }y t } 2 ˘1{2
, where p a n,t and x M n,t denote the estimated abundances and EMs. To evaluate the EMs, we also computed the average spectral angle mapper (SAM) as SAM M " 1

T N P ř T t"1 ř N n"1 ř P j"1 arccos
´mJ n,t,j x mn,t,j }mn,t,j }}x mn,t,j } ¯, in which m n,t,j and x m n,t,j are the true and estimated EM signatures for time t, pixel n and EM j.

A. Simulations with synthetic data

Two synthetic datasets were considered with spatiotemporal abundance and EM variability. The first dataset, referred to as Sequence 1 (DS1), contains P " 3 EMs and T " 6 HIs. The HIs are generated from sequences of abundance maps with N " 50 ˆ50 pixels containing localized abrupt changes for t P t2, 3, 4, 5u (depicted in the first row of Figure 5). The EMs for each pixel and time instant were generated as follows. First, three signatures with L " 224 bands were selected from the USGS library and used as a reference EM matrix. Then, for the first time instant (t " 1), spatial EM variability was introduced by following the model in [START_REF] Thouvenin | Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[END_REF], in which the EMs in each pixel (M n,1 , n " 1, . . . , N ) were generated by multiplying the reference signatures with piecewise linear random scaling factors with amplitude in the interval r0.85, 1.15s. For each subsequent time instant t ą 1, the EMs were also generated as scaled versions of the reference spectral signatures. However, to introduce temporal EM variability, the scaling factors at time t are defined to be the sum of the scaling factors at time t ´1 plus random piecewise linear functions in the range r´0.1, 0.1s. Samples of the generated EMs can be seen in Figure 6. These EM matrices M n,t are then used to generate the HI pixels using the LMM (1), with the measurement noise r n,t being white and Gaussian with an SNR of 30 dB. The second dataset, referred to as Data Sequence 2 (DS2), contained P " 4 EMs and N " 50 ˆ50 pixels. A sequence of abundance maps generated randomly according to a Gaussian random field and containing small, spatially compact abrupt changes was considered to generate T " 15 HIs. To introduce realistic spectral variability, the EM signatures at each pixel and time instant were randomly selected from a set of pure pixels of water, vegetation, soil and road that were manually extracted from the Jasper Ridge HI, with L " 198 bands. The HI sequence was then generated according to the multitemporal LMM (1), with the r n,t being white Gaussian noise with an SNR of 30 dB. The parameters of the ReDSUNN were K " 10 and σ ψ " 10 ´5 for DS1, and K " 2 and σ ψ " 10 ´5 for DS2. The quantitative results are presented in Table II, while the visual results (only shown for DS1 due to space limitations) are depicted in Figures 5 and7.

1) Discussion: It can be seen from results, while the performance of the Kalman filter was good for DS1 but very poor for DS2. The FLCS, which does not take temporal information of spatial EM variability into account, did not perform very well, having the worse abundance reconstructions on average for both datasets. From the estimated abundances in Figure 5, it can be seen that ReDSUNN's results are the closest to the ground truth. However, the results for all methods were relatively noisy. The abundances recovered by the Kalman filter, OU and HBUN indicated more heavily mixed pixels. FCLS achieves reasonable performance for t ď 4, but led to a completely wrong estimation for t " 6.

The changes occurring in the ground truth abundances can be observed in the estimations of all methods, although they are more clearly visible in the Kalman and ReDSUNN results since these methods led to a larger separation between the different materials. The visual abundance results for DS2 (not shown due to space limitations) were qualitatively similar to those of DS1, with the exception that the Kalman filter failed to identify the soil EM for all images in the sequence, which explains its poor performance.

The ReDSUNN method also obtained the best EM estimation performance all metrics except for the SAM in DS2, in which OU achieved the best result followed by HBUN. The Kalman filter obtained good results for DS1 (close to ReDSUNN), but poor results in DS2. This happened despite the Kalman filter obtaining reasonable image reconstruction errors NRMSE Y for both DS1 and DS2. Samples of the true and estimated EMs in Figures 6 and7 (only shown for DS1 and for ReDSUNN due to space limitations) indicate that the EMs are correctly recovered. However, there are some differences, particularly in the shape of the first EM (which show higher amplitude in the ground truth compared to the estimates). Moreover, the amount of variability was lower in The lowest image reconstruction errors (NRMSE Y ) were obtained by OU and HBUN, while those obtained by ReDSUNN were similar to those by FCLS. This is expected, since NRMSE Y is closely related to the number of learnable parameters of each algorithm, and is not directly related to the abundance or EM reconstruction performance. This explains the higher reconstruction error by ReDSUNN since, as discussed in Section V-D, the shared parametrization of the variational posterior PDF leads to a relatively low number of parameters, which also helps to mitigate overfitting. Nonetheless, for the synthetic data sequences (DS1 and DS2) ReDSUNN still has between 30% and 50% more learnable parameters than OU. Its parametrizaition becomes significantly more favorable when the images have a larger amounts of pixels, such as in the experiments with the Lake Tahoe images presented in Section VI-C. The computation times show a clear separation between FCLS and OU, which were faster, and HBUN, the Kalman filter and ReDSUNN, which took longer to run. This indicates that the proposed method has a competitive computational performance when compared to more complex algorithms.

B. Sensitivity

To measure the influence of different hyperparameters on the performance of the method, we evaluated how NRMSE A varied as a function of the hyperparameters, namely, the number of basis vectors for the variability model, K, and the innovation standard deviation of the EM variability parameters, σ ψ . The results for both DS1 and DS2 can be seen in Figure 8. It can be seen that for DS1, the performance of ReDSUNN is not heavily affected by the number of bases K within the evaluated range. However, σ ψ has a larger impact on the result, with smaller values leading to a lower NRMSE A . For DS2, smaller values for both parameters generally lead to lower NRMSE A results, although there performance varied more with K and σ ψ . For both datasets, small variations of these parameters around the optimal values lead to similar results. In general, the larger the value of σ ψ , the more temporal EM variability is allowed by the model, whereas the larger the value of K, the more complex the spatial and temporal EM variability the model can represent. Devising a methodology to automatically tune these parameters is an interesting question for future work.

C. Simulations with real data

To evaluate the performance of the algorithms on real data, we considered the Lake Tahoe HI sequence, which was originally described in [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF]. It consists of a sequence of T " 6 images acquired over the Lake Tahoe area by the AVIRIS instruments, which are depicted in true color in Figure 9. Each HI contained N " 16500 pixels, and L " 173 bands were left after the removal of low-SNR and water absorption bands. This scene contains P " 3 predominant EMs, consisting of soil, water and vegetation, and considerable changes on the lake and on the crop circles can be observed between the images. The parameters of the ReDSUNN were set as K " 3, σ ψ " 1, while the parameters of OU, HBUN and of the Kalman filter were selected as described in their original publications. The recovered abundances are depicted in Figure 11, while the recovered EMs (only shown for ReDSUNN due to space limitations) are shown in Figure 10. The reconstruction errors and the processing times are presented in Table III.

1) Discussion: From Figure 11, it can be seen that the FCLS method did not achieve a good performance in general, particularly for the fifth image in which there was a considerable confusion between the soil and vegetation EMs. The remaining algorithms achieve more stable performance due to taking the temporal information into account. The OU and the HBUN algorithms (both of which use the PLMM [START_REF] Thouvenin | Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[END_REF] model to represent the temporal endmember variability) behaved similarly to each other. Although these methods performed more stably than the FCLS, they still presented considerable water abundances outside of the lake region, which are predominantly composed by soil. The performance of the Kalman filter method was relatively poor for the third, fourth and sixth images (in which the area of the lake is small), and contained a considerable amount of artifacts. This happens since the Kalman filter method assumes the abundances to be constant over time when estimating the EMs. Consequently, it is not able to handle large abundance changes in the image sequence. The abundances estimated by ReDSUNN, on the other hand, showed a clear separation between the different materials, adequately capturing the abundance changes occurring in the HIs. Moreover, larger concentrations of mixed pixels were observed in regions that are meaningful, such as at the drying edge of the lake in the third image, and in some parts of the crop circles. The EMs recovered by ReDSUNN show considerable variability in soil and vegetation spectral, particularly over space, while the water spectra shows little variability. Moreover, the variability of the EMs can be spectrally localized, which can be observed most clearly in the temporal signatures of soil. Moreover, spatial EM variability was more significant than temporal EM variability. Note that the EM variability was also more significant in this example compared to the experiments with synthetic data since a larger value for the hyperparameter σ ψ was selected.

The results in Table III show that HBUN and OU obtain the smallest reconstruction errors (NRMSE Y ), while those obtained by Kalman and ReDSUNN, which have less learnable parameters, were larger, with FCLS having the largest NRMSE Y . The ratio between the computation times were similar to the synthetic examples, with the Kalman filter being the slowest and the OU the fastest among the MTHU methods that account for temporal information, and the proposed ReDSUNN method achieving intermediate results. This indicates that the methods scale similarly with the image size. Nevertheless, developing more efficient MTHU algorithms is an interesting subject for future work.

VII. CONCLUSIONS

This paper proposed a multitemporal hyperspectral unmixing method based on a variational recurrent neural network. A low-dimensional, dynamical state space model was presented to represent the spatial and temporal variations of the endmember spectra by expanding it over a small set of spectrally smooth basis vectors. The dynamics of the abundances were modelled using a Dirichlet distribution, which was approximated as a Gaussian in the softmax basis in order to improve the efficiency of the inference process. Based on this generative model, variational inference was considered to perform unmixing by approximating the posterior distribution of the abundances and endmembers. The Markov and independence properties of the model were also used to improve the efficiency of the solution. The posterior distribution was parameterized using a combination of a simple, physically interpretable, model and LSTM recurrent neural networks to improve flexibility while maintaining the physical interpretability of the abundances. In the proposed framework, all parameters were computed using stochastic backpropagation. Experimental results indicate that the proposed algorithm achieves better unmixing performance when compared to state-of-the-art methods, at a similar computational complexity, using both synthetic and real datasets. [START_REF] Uezato | Hyperspectral unmixing with spectral variability using adaptive bundles and double sparsity[END_REF] Due to the (conditionally) Gaussian assumptions in the generative model and in the variational posterior, the three terms inside the expectations in [START_REF] Uezato | Hyperspectral unmixing with spectral variability using adaptive bundles and double sparsity[END_REF] can be computed analytically. The first term in [START_REF] Uezato | Hyperspectral unmixing with spectral variability using adaptive bundles and double sparsity[END_REF] is the log-likelihood of a Gaussian PDF, which can be computed from [START_REF] Imbiriba | Hybrid Neural Network Augmented Physics-based Models for Nonlinear Filtering[END_REF]. The second and third terms are KL divergences between Gaussian PDFs, which can be computed using the general result for two Gaussians of dimension D, given by KL `N pµ 1 , Σ 1 q › › N pµ 2 , Σ 2 q ˘" 1 2 `log |Σ2| |Σ1| ´D `trtΣ ´1 2 Σ 1 u `pµ 2 ´µ1 q J Σ ´1 2 pµ 2 ´µ1 q ˘. The second (resp., third) term are thus computed substituting the mean and covariance from [START_REF]Real-time spatiotemporal spectral unmixing of MODIS images[END_REF] and ( 20), (21) (resp., [START_REF] Dudley | A multi-temporal spectral library approach for mapping vegetation species across spatial and temporal phenological gradients[END_REF] and ( 18), [START_REF]Multi-temporal hyperspectral mixture analysis and feature selection invasive species mapping in rainforests[END_REF]) in the expression above by using the fact that c n,t and ψ n,t are independent in the generative model. For more details, see, e.g., [START_REF] Doersch | Tutorial on variational autoencoders[END_REF].

APPENDIX A COMPUTING THE TERMS IN

Fig. 3 :

 3 Fig. 3: Illustrative diagram of the proposed generative model.

Fig. 5 :

 5 Fig. 5: Estimated abundance maps and ground truth for DS1, shown as composite color maps (that is, the abundances of the EMs #1, #2 and #3 are represented as the red, green, and blue color channels).

Fig. 6 :Fig. 7 :

 67 Fig. 6: True EMs for the DS1, sampled over space, for time instant t " 3 (top), and over time, for pixel n " 1 (bottom). EM #1 EM #2 EM #3

Fig. 8 :

 8 Fig. 8: Abundance RMSE as a function of hyperparameters K and σ ψ for DS1 (left) and DS2 (right).

  the retrieved EMs compared to the ground truth; this occurs for the synthetic examples since the hyperparameter σ ψ , which controls the flexibility of the EM model, was selected to provide the best performance in terms of NRMSE A , leading to relatively small σ ψ values. This is further illustrated in the experiments shown in Section VI-B.

Fig. 9 :Fig. 10 :

 910 Fig. 9: True color depiction of the Lake Tahoe HIs and their acquisition dates. Water Soil Vegetation

Fig. 11 :

 11 Fig. 11: Estimated abundances for the water (left panel), soil (central panel) and vegetation (right panel) EMs of the Lake Tahoe HIs.

  Experimental results with synthetic and real data show that ReDSUNN outperforms state of the art MTHU algorithms. Codes are available at https://github.com/ricardoborsoi/ReDSUNN.

formulation of model-based MTHU with both spatial and temporal EM variability, solved using stochastic backpropagation; ' a parametrization of the posterior distribution of the abundances and EMs combining bidirectional RNNs and a physically interpretable model. The proposed method is called ReDSUNN for Recurrent hyperSpectral Unmixing with Neural Networks.

  ) ; Compute p cn,t, p ψ n,t as the means of q φ pct, ψ t |p c t´1 , p ψ t´1 , yq using (44) and (45) ; M n,t, for n " 1, . . . , N and t " 1, . . . , T .

	6	for t " 1, . . . , T do
	7
	8	end
	9 end
	10 Set p an,t " πpp cn,tq, x M n,t " x M 0 d p1 `D vec ´1p p ψ n,t qq ;
	Output: p an,t, x

TABLE I :

 I Variables to be estimated and number of parameters.

	Generative model (θ)

TABLE II :

 II Quantitative results of the simulations using synthetic data.NRMSE A NRMSE M SAM M NRMSE Y Time

			Data Sequence 1 -DS1
	FCLS	0.537	-	-	0.086	2.7
	OU	0.434	0.342	0.260	0.051	24.9
	HBUN	0.479	0.355	0.162	0.050	542.6
	Kalman	0.356	0.124	0.076	0.061	2422.8
	ReDSUNN	0.318	0.117	0.075	0.089	479.0
			Data Sequence 2 -DS2
	FCLS	0.500	-	-	0.122	7.3
	OU	0.335	0.256	0.120	0.055	60.6
	HBUN	0.474	0.515	0.141	0.050	2166.0
	Kalman	0.659	12.222	0.496	0.108	5937.4
	ReDSUNN	0.294	0.203	0.289	0.160	1231.3

TABLE III :

 III Quantitative results for the Lake Tahoe HI sequence.

		FCLS	OU	HBUN Kalman ReDSUNN
	NRMSE Y	0.321	0.058	0.054	0.185	0.114
	Time	16.1	92.5	3381.9	4607.7	2857.2

This work has been partially supported by the NSF under Award ECCS-1845833.