Generative Simplex Mapping: Non-Linear Endmember Extraction and Spectral Unmixing for Hyperspectral Imagery

Abstract

We introduce a new model for non-linear endmember extraction and spectral unmixing of hyperspectral imagery called Generative Simplex Mapping (GSM). The model represents endmember mixing using a latent space of points sampled within a $(n-1)$-simplex corresponding to n unique sources. Barycentric coordinates within this simplex are naturally interpreted as relative endmember abundances satisfying both the abundance sum-to-one and abundance non-negativity constraints. Points in this latent space are mapped to reflectance spectra via a flexible function combining linear and non-linear mixing. Due to the probabilistic formulation of the GSM, spectral variability is also estimated by a precision parameter describing the distribution of observed spectra. Model parameters are determined using a generalized expectation-maximization algorithm, which guarantees non-negativity for extracted endmembers. We first compare the GSM against three varieties of non-negative matrix factorization (NMF) on a synthetic data set of linearly mixed spectra from the USGS spectral database. Here, the GSM performed favorably for both endmember accuracy and abundance estimation with all non-linear contributions driven to zero by the fitting procedure. In a second experiment, we apply the GTM to model non-linear mixing in real hyperspectral imagery captured over a pond in North Texas. The model accurately identified spectral signatures corresponding to near-shore algae, water, and rhodamine tracer dye introduced into the pond to simulate water contamination by a localized source. Abundance maps generated using the GSM accurately track the evolution of the dye plume as it mixes into the surrounding water.

1. Introduction

Hyperspectral imaging has emerged as a keystone technology in remote sensing, where the ability to discern variations between high-resolution spectra supports a plethora of critical applications such as environmental monitoring, biodiversity conservation, sustainable agriculture, and more. In recent years, many remote sensing platforms have been deployed with hyperspectral imaging payloads such as the Italian PRISMA mission launched in 2019, the German EnMAP launched in 2022, and recently, NASA’s PACE satellite launched in 2024 [1,2,3]. Many future missions also plan to incorporate hyperspectral imaging capabilities, such as the European Space Agency’s CHIME, which will include over 200 bands spanning visible, near-infrared (NIR), and short-wave infrared (SWIR) wavelengths [4]. The continued development of hyperspectral imaging technology has also led to a considerable reduction in size that enables its inclusion in the payloads of small unmanned aerial vehicles (UAVs) [5,6]. Despite the proliferation of hyperspectral imaging data sources, the considerable increase in data volume associated with hyperspectral images (HSIs) poses significant challenges to real-time analysis at scale.

Many approaches have been developed to make sense of the HSI data. For example, spectral indices such as the popular normalized difference vegetation index (NDVI) can be computed by taking ratios of spectral bands tailored to track specific reflectance characteristics [7,8]. These indices have the advantage of being easy to compute, but suffer from significant variability between instruments while ignoring most of the information captured in HSI spectra [9]. An alternative approach is to pair HSI data with in situ measurements to enable supervised models that map spectra directly to parameters of interest. However, this approach relies on serendipitous satellite overpasses above sensing sites to generate sufficient quantities of aligned data for model training and evaluation. For example, Aurin et al. combined data from over 30 years of oceanographic field campaigns with paired satellite imagery to develop robust models for the inversion of key water quality indicators such as colored dissolved organic matter [10]. This approach can be accelerated by combining UAV-based hyperspectral imaging with rapid in situ data collection using autonomous boats [11,12]. However, these supervised methods rely on a priori knowledge of expected sources to identify appropriate reference sensors for data collection. In light of these limitations, unsupervised methods are needed that can reliably extract source signatures from high-dimensional HSIs.

The spatial resolution of hyperspectral images generally results in pixels with mixed signals from multiple sources, called endmembers. The task of unsupervised source identification using HSI data therefore involves two steps: endmember extraction and abundance estimation. Techniques such as vertex component analysis (VCA), the pixel purity index (PPI), and N-FINDR solve this first task by identifying HSI endmember spectra assuming the presence of some pure (unmixed) pixels [13,14,15]. By further assuming a linear mixing model (LMM) in which observed spectra are described by a linear combination of endmembers with non-negative abundances, HSIs can be unmixed using a variety of techniques such as constrained least squares and multi-objective optimization [16,17,18]. Among these methods, non-negative matrix factorization (NMF) is a widely used approach that extracts endmember spectra and unmixes abundances simultaneously via matrix factorization [19,20]. The update equations for NMF can be formulated as multiplicative updates, which guarantee the non-negativity of endmember spectra and their associated abundances [21]. For this reason, the continued development of new NMF varieties remains an active area of research.

In realistic scenes, multiple scattering and surface variability can easily challenge the assumption of linear mixing [22]. Water-based HSIs specifically are prone to non-linear mixing effects due to absorption features of dissolved and suspended substances, fluorescence of organic matter, and particulate scattering in turbid waters [23,24,25]. With the growing popularity of deep learning approaches in remote sensing, a variety of models based on autoencoder architectures have been introduced for unmixing HSI data [26,27,28,29]. However, the complexity introduced by these models significantly impacts training time and decreases the interpretability of the resulting model. An ideal approach should enable both endmember extraction and non-linear unmixing while accounting for spectral variability.

The self-organizing map (SOM) is an unsupervised machine learning method that maps high-dimensional data to a low-dimensional grid while preserving the topological relationships between data points [30]. This low-dimensional representation provides a convenient way to visualize HSI data, while the weight vectors for each SOM node can be interpreted as representative spectra [31,32,33]. If labeled reference spectra are available, the SOM can be used to enable semi-supervised labeling of HSI pixels [34]. The SOM has also been shown to be effective for the compression of HSIs acquired by a CubeSat [35]. Despite these capabilities, the SOM does not offer a probabilistic interpretation and relies on a heuristic training procedure with hyperparameters that can be challenging to tune. To address these shortcomings, Bishop et al. introduced Generative Topographic Mapping (GTM), which is a probabilistic latent-variable model inspired by the SOM [36]. When choosing a two-dimensional latent space, the GTM can be used to visualize the distribution of HSI spectra, while the non-linear mapping from latent space nodes to the HSI data space provides endmembers [37]. Unfortunately, the rectangular latent space grid employed by the GTM does not directly translate into endmember abundances, and the expectation-maximization (EM) algorithm used to train the GTM does not guarantee non-negativity of GTM node spectra.

In this paper, we introduce a new variant of the GTM, dubbed Generative Simplex Mapping (GSM), which can extract endmember spectra and unmix non-linear mixtures. By replacing the rectangular latent space of the GTM with a gridded $(n-1)$-simplex, the vertices of the GSM can be immediately interpreted as endmembers corresponding to n unique spectral signatures. The mapping from the latent space to the HSI data space accounts for linear and non-linear spectral mixing while barycentric coordinates for the latent space simplex represent relative endmember abundances. Additionally, by taking inspiration from the multiplicative updates of NMF, the GSM algorithm maintains the non-negativity of resulting endmember spectra. If only linear mixing is present, the GSM algorithm drives non-linear contributions to 0. Prior distributions included for GSM model weights yield hyperparameters, which can be tuned to control the smoothness of the resulting spectra and the degree of non-linear mixing applied. In summary, the key innovations introduced by the GSM are as follows:

• The GSM can model linear and non-linear spectral mixing;
• The GSM does not assume the presence of pure pixels in the data set;
• The probabilistic formulation of the GSM accounts for spectral variability;
• The simplex used for the latent space structure of the GSM is directly interpretable and forces abundances to satisfy both the abundance sum-to-one and abundance non-negativity constraints;
• The fitting procedure introduced for the GSM maintains non-negativity of endmember spectra.

The remainder of this paper is structured as follows. Section 2 describes the proposed method. Section 3 and Section 4 describe experiments using simulated and real HSI data to evaluate the GSM model. Section 5 discusses additional applications and extensions of the GSM to be explored in future work. Finally, Section 6 finishes the paper with some closing remarks.

2. Generative Simplex Mapping

In the original GTM formulation, the data vectors $\mathbf{x}$ (reflectance spectra of length D) are described by latent variables $\mathbf{z}$ mapped into the data space by a non-linear function $\psi$. The data space distribution is taken to be normal with a precision parameter $\beta$ to account for measurement noise and spectral variability. The GSM uses same structure, that is

$$p(\mathbf{x}\mid \mathbf{z},\mathbf{W},\beta )={\left({\displaystyle \frac{\beta }{2\pi }}\right)}^{D/2}exp\left(-{\displaystyle \frac{\beta }{2}}{\parallel \psi (\mathbf{z};\mathbf{W})-\mathbf{x}\parallel }^2\right)$$

(1)

where $\mathbf{z}=(z_1,z_2,...,z_{N_v})$ corresponding to $N_v$-many sources and $\mathbf{W}$ is a $D\times M$ matrix of model weights that parameterize the mapping $\psi$. As Equation (1) indicates, the precision parameter corresponds to a standard deviation of $\sqrt{{\beta }^{-1}}$.

Assuming data are uniformly sampled from an embedded manifold in the data space, the GTM models the latent space using a rectangular grid with K-many nodes of equal prior probability. To adapt the GTM to describe endmember mixing, the GSM makes two key changes. The first is to replace the rectangular GTM grid with a gridded simplex with $N_v$ vertices. Barycentric coordinates for each node then describe the relative abundance of each endmember with simplex vertices corresponding to pure sources. The second change is to augment the equal prior probabilities of the GTM latent space with adaptive mixing coefficients ${\pi }_k$ to allow for nonuniform sampling across the abundance simplex. This addition allows the GSM to flexibly address many possible mixing scenarios without an a priori assumption for a specific mixing distribution. Together, these changes result in a prior distribution for the latent space given by

$$p\left(\mathbf{z}\right)=\sum _k^K{\pi }_k\delta (\mathbf{z}-{\mathbf{z}}_k)\mathrm{where}\sum _k^K{\pi }_k=1.$$

(2)

where $\delta (·)$ is the Dirac delta function and ${\mathbf{z}}_k$ are the positions of each node within the simplex.

For a data set containing N-many records, these definitions yield a log likelihood function given by

$$\mathcal{L}=\sum _n^{N}ln\left(\sum _k^K{\pi }_{k}p({\mathbf{x}}_n\mid {\mathbf{z}}_k,\mathbf{W},\beta )\right)$$

(3)

which can be maximized to obtain optimal values for model weights $\mathbf{W}$, mixing coefficients ${\pi }_k$, and precision $\beta$. Rather than optimizing Equation (3) directly, we instead choose a particular form for $\psi$ to allow fitting the GSM via an EM algorithm.

In the standard LMM model, the mapping $\psi$ is given by $\psi (\mathbf{z};\mathbf{W})=\mathbf{W}\mathbf{z}$ where the columns of $\mathbf{W}$ correspond to endmember spectra. To model non-linear mixing, $\mathbf{z}$ is replaced by the output of M-many activation functions such that $\psi (\mathbf{z};\mathbf{W})=\mathbf{W}\varphi \left(\mathbf{z}\right)$. The activations applied to each GSM node can then be collected to form a matrix with elements ${\mathrm{\Phi }}_{km}={\varphi }_m\left({\mathbf{z}}_k\right)$. For $m\le N_v$, we take ${\mathrm{\Phi }}_{km}={\left[{\mathbf{z}}_k\right]}_m$ to model linear mixing. The remaining $M-N_v$ activations are computed using radial basis functions (RBFs) with centers ${\mu }_m$ distributed throughout the simplex (but not at the vertices) with

$${\mathrm{\Phi }}_{km}=\left\{\begin{array}{cc}{\displaystyle \frac{s-\parallel {\mathbf{z}}_k-{\mu }_m\parallel }{s}};\hfill & \parallel {\mathbf{z}}_k-{\mu }_m\parallel \le s\hfill \\ 0;\hfill & \parallel {\mathbf{z}}_k-{\mu }_m\parallel >s\hfill \end{array}\right.$$

(4)

where s is the spacing between RBF centers. A visual representation of the non-linear activation functions in Equation (4) is shown below in Figure 1 for a 2-component GSM.

In this form, the first $N_v$ columns of $\mathbf{W}$ correspond to endmember spectra, while the remaining columns account for additional non-linear effects. Equation (4) is specifically chosen so that no non-linear contributions are possible for pure spectra at the vertices of the simplex. At all other points, the output of $\psi$ involves both linear and non-linear contributions. If only linear mixing is present, the GSM training algorithm should therefore drive $W_{dm}$ to 0 for $m\ge N_v$. A visualization of the GSM model for three sources is shown in Figure 2.

To further constrain the model, we introduce prior distributions on the weights $\mathbf{W}$. For $m\le N_v$ we take $W_{dm}\sim \mathcal{N}(0,{\lambda }_e^{-1})$, corresponding to a zero-mean Gaussian with variance ${\lambda }_e^{-1}$. For $m>N_v$ we use a zero-mean Laplace distribution, $W_{dm}\sim {\displaystyle \frac{{\lambda }_w}{2}}exp(-{\lambda }_w\left|W_{dm}\right|)$, with scale parameter ${\lambda }_w^{-1}$. Under these choices, ${\lambda }_e$ corresponds to $L_2$ regularization on endmember spectra while ${\lambda }_w$ corresponds to $L_1$ regularization on the non-linear activations. In other words, ${\lambda }_e$ governs the smoothness of the resulting endmembers, while ${\lambda }_w$ encourages sparsity for the non-linear contributions.

An EM algorithm for the GSM model can now be formulated as follows. Suppose that we have current estimates for the model weights $\mathbf{W}$, mixing coefficients ${\pi }_k$, and precision parameter $\beta$. During the expectation step we compute the posterior probabilities, that is, the responsibility of each GSM node for each spectrum in the data set:

$$R_{kn}=p({\mathbf{z}}_k\mid {\mathbf{x}}_n,\mathbf{W},\beta )={\displaystyle \frac{{\pi }_{k}p({\mathbf{x}}_n\mid {\mathbf{z}}_k,\mathbf{W},\beta )}{{\displaystyle \sum _{k^{\prime }}^K}{\pi }_{k^{\prime }}p({\mathbf{x}}_n\mid {\mathbf{z}}_{k^{\prime }},\mathbf{W},\beta )}}.$$

(5)

For the maximization step, we consider the expectation of the penalized complete-data log likelihood given by

$$\begin{array}{cc}\hfill Q& =\sum _n^N\sum _k^K{R}_{kn}\left(ln{\pi }_k+{\displaystyle \frac{D}{2}}ln\left({\displaystyle \frac{\beta }{2\pi }}\right)-{\displaystyle \frac{\beta }{2}}\sum _d^D{\left(\sum _m^M{W}_{dm}{\mathrm{\Phi }}_{km}-X_{nd}\right)}^2\right)\hfill \\ \hfill & +{\displaystyle \frac{N_{v}D}{2}}ln\left({\displaystyle \frac{{\lambda }_e}{2\pi }}\right)-{\displaystyle \frac{{\lambda }_e}{2}}\sum _d^D\sum _{m=1}^{N_v}{W}_{dm}^2\hfill \\ \hfill & +(M-N_v)Dln\left({\displaystyle \frac{{\lambda }_w}{2}}\right)-{\lambda }_w\sum _d^D\sum _{m=N_v+1}^M{W}_{dm}\hfill \end{array}$$

(6)

where $X_{nd}$ is the d-th component of the n-th spectrum in the data set. Equation (6) is then maximized with respect to ${\pi }_k$, $\beta$, and $\mathbf{W}$ to obtain new parameter values. For a detailed overview of the EM procedure, the reader is directed to ref. [38].

For ${\pi }_k$, optimization can be performed using Lagrange multipliers to maintain the condition that ${\sum }_k{\pi }_k=1$. Doing so yields

$${\pi }_k^{\mathrm{new}}={\displaystyle \frac{1}{N}}\sum _n{R}_{kn}$$

(7)

Optimization of Equation (6) with respect to $\mathbf{W}$ leads to a linear system that can be solved using standard numerical methods. However, in this form, we cannot guarantee the non-negativity of $\mathbf{W}$ required to describe reflectance spectra. Therefore, we take inspiration from the multiplicative updates for NMF introduced by Lee and Seung [19]. A standard gradient-based update for $\mathbf{W}$ would normally take the form

$${\mathbf{W}}_{\mathrm{new}}=\mathbf{W}+\eta {\displaystyle \frac{\partial Q}{\partial \mathbf{W}}}$$

(8)

for some learning rate $\eta$. Therefore, we differentiate to obtain

$${\displaystyle \frac{\partial Q}{\partial \mathbf{W}}}=-\beta \mathbf{W}{\mathbf{\Phi }}^T\mathbf{G}\mathbf{\Phi }-\mathrm{\Lambda }+\beta {\mathbf{X}}^T{\mathbf{R}}^T\mathbf{\Phi }$$

(9)

where $\mathbf{G}$ is a diagonal matrix with $G_{kk}={\sum }_n{R}_{kn}$ and $\mathrm{\Lambda }$ is given by

$${\mathrm{\Lambda }}_{dm}=\left\{\begin{array}{cc}{\lambda }_e{W}_{dm};\hfill & m\le N_v\hfill \\ {\lambda }_w;\hfill & m>N_v\hfill \end{array}\right.$$

(10)

If we allow individual learning rates ${\eta }_{dm}$ for each element of $\mathbf{W}$, then choosing

$${\eta }_{dm}={\displaystyle \frac{W_{dm}}{{\left(\beta \mathbf{W}{\mathbf{\Phi }}^T\mathbf{G}\mathbf{\Phi }\right)}_{dm}+{\mathrm{\Lambda }}_{dm}}}$$

(11)

results in a multiplicative update rule given by

$$W_{dm}^{new}=W_{dm}·{\displaystyle \frac{{\left(\beta {\mathbf{X}}^T{\mathbf{R}}^T\mathbf{\Phi }\right)}_{dm}}{{\left(\beta \mathbf{W}{\mathbf{\Phi }}^T\mathbf{G}\mathbf{\Phi }\right)}_{dm}+{\mathrm{\Lambda }}_{dm}}}$$

(12)

From Equation (12), it is clear that we are multiplying $W_{dm}$ by strictly non-negative values, and therefore, non-negative $W_{dm}$ will remain so during each update. This update can also be repeated multiple times during each M-step to accelerate convergence.

Optimizing Equation (6) with respect to $\beta$ yields the final update equation:

$${\displaystyle \frac{1}{{\beta }^{\mathrm{new}}}}={\displaystyle \frac{1}{ND}}\sum _n^N\sum _k^K{R}_{kn}{\parallel \psi ({\mathbf{z}}_k;\mathbf{W})-{\mathbf{x}}_n\parallel }^2.$$

(13)

To train a GSM model, weights $\mathbf{W}$ are randomly initialized to positive values. The mixing coefficients are initially set to ${\pi }_k=1/K$. Finally, the precision parameter $\beta$ is initialized to the variance of the $(N_v+1)$-th principal component. The matrix $\mathbf{\Phi }$ is pre-computed and stored for reuse during each EM cycle together with the data matrix $\mathbf{X}$. After initialization, the expectation and maximization steps are repeated in turn until Q converges to a predetermined tolerance level. For large $N_v$, we note that generating a regular grid within the simplex becomes cumbersome as the number of grid nodes scales as $\left({\displaystyle \genfrac{}{}{0pt}{}{k+N_v-2}{N_v-1}}\right)$ for k nodes per edge leading to a large responsibility matrix. An alternative approach is to randomly sample points within the simplex to obtain a total of K nodes. A Dirichlet distribution

$$p\left(\mathbf{z}\right)={\displaystyle \frac{\mathrm{\Gamma }\left({\sum }_i^{N_v}{\alpha }_i\right)}{{\prod }_i^{N_v}\mathrm{\Gamma }\left({\alpha }_i\right)}}\prod _i^{N_v}{\left[\mathbf{z}\right]}_i^{{\alpha }_i-1}$$

(14)

with all ${\alpha }_i=1$ can be used to uniformly sample within the simplex. Since the mixing coefficients ${\pi }_k$ are adaptive, variability in node separation should not significantly impact the resulting GSM.

The probabilistic form of the GSM means that a variety of information criteria can be used to evaluate the model fits. In this paper we consider two metrics, the Bayesian Information Criterion (BIC),

$$\mathrm{BIC}=Pln\left(N\right)-2\mathcal{L},$$

(15)

and the Akaike Information Criterion (AIC),

$$\mathrm{AIC}=2P-2\mathcal{L}$$

(16)

where P is the total number of model parameters and $\mathcal{L}$ is the log likelihood from Equation (3).

The map $\psi$ provides the representation of each node $z_k$ in the data space. Importantly, applying $\psi$ to each vertex extracts endmember spectra from the GSM. Slices $R_{\left[:,n\right]}$ of the matrix $\mathbf{R}$ define the responsibility of each latent node $z_k$ for the n-th spectrum $x_n$ in the data set. Therefore, once the GSM has been trained, $\mathbf{R}$ can be used to unmix endmember abundances by representing each record in the latent space via the mean:

$${\widehat{\mathbf{z}}}_n=\sum _k^K{R}_{kn}{\mathbf{z}}_k.$$

(17)

A freely available implementation of the GSM is provided in [39]. The code is written in the Julia programming language and follows the Machine Learning in Julia (MLJ) common interface [40,41].

3. Experiments

3.1. Linear Mixing: Comparison to NMF

To illustrate the effectiveness of the GSM, we first demonstrate its ability to model linear mixing. This serves as an important limiting case, since linearly mixed spectra should not lead to the spurious introduction of non-linear contributions by the GSM. The goal of this first test is therefore to demonstrate that the GSM drives non-linear weights to zero for linearly mixed data while providing a fair test to compare the GSM to a well-established linear mixing model. This ability clearly distinguishes the GSM from other non-nonlinear unmixing approaches such as autoencoders, which, by their design, include non-linear mixing even when it is not present in the underlying data. To this end, a synthetic data set comprising linear mixtures of three sample spectra from the U.S. Geological Survey digital spectral library was generated by sampling 1000 abundance vectors from a Dirichlet distribution with ${\alpha }_1={\alpha }_2={\alpha }_3=1/3$ [42]. These source spectra and abundances are visualized in Figure 3.

Zero-mean Gaussian noise was added to the data to yield 9 data sets with signal-to-noise ratios (SNRs) ranging from 0 to ∞ to examine the impact of random noise on GSM performance. For each of these data sets, a GSM was trained using 25 nodes per edge with ${\lambda }_e$ set to the default value of $0.01$ and ${\lambda }_w$ fixed to 100. The large value of ${\lambda }_w$ was chosen to encourage the GSM to prefer models with limited non-linear mixing.

For comparison, three NMF models were also trained on each data set, as this method is both highly popular and does not include the pure pixel assumption common to other techniques like VCA and PPI. Countless variations on the original NMF method have been introduced into the literature, with one review identifying more than 100 distinct NMF varieties [20]. For the purpose of comparing to the GSM, we considered the standard ${\ell }_2$ and KL-divergence formulations introduced by Lee and Seung [21]. We also included the robust ${\ell }_{2,1}$ NMF as described by Kong et al. [43]. We note that the goal of this test is not to prove the GSM is superior to other models for linear mixing, but rather to demonstrate that the GSM can model linearly mixed data without introducing unnecessary complexity.

Four metrics were used to compare model performance. For endmember extraction, the mean spectral angle and mean RMSE between true endmembers ${\rho }_i$ and extracted endmembers ${\widehat{\rho }}_i$ were computed where the spectral angle for the i-th endmember is defined as

$$\theta ({\rho }_i,{\widehat{\rho }}_i)=arccos\left({\displaystyle \frac{\langle {\rho }_i,{\widehat{\rho }}_i\rangle }{\parallel {\rho }_i\parallel ·\parallel {\widehat{\rho }}_i\parallel }}\right)$$

(18)

and the RMSE for the i-th endmember is

$$\mathrm{RMSE}({\rho }_i,{\widehat{\rho }}_i)=\sqrt{{\displaystyle \frac{1}{D-1}}\sum _d^D{\left({\rho }_i\left({\lambda }_d\right)-{\widehat{\rho }}_i\left({\lambda }_d\right)\right)}^2}.$$

(19)

Abundance estimation was similarly evaluated using the mean RMSE between true abundance and estimated abundance values for each endmember. Finally, the reconstruction RMSE, that is, the RMSE computed between the original data set and its reconstruction via the extracted endmembers and abundances, was computed. This provides a model-agnostic criterion to guarantee that each model sufficiently converged during training.

3.2. Non-Linear Mixing: Water Contaminant Identification

To assess the ability of the GSM to unmix realistic scenes likely to involve non-linear mixing effects, we consider a data set of real HSI collected in Montague, North Texas on 9 December 2020. A Freefly Alta-X autonomous quadcopter was used as a UAV platform and equipped with a Resonon Pika XC2 visible+near-infrared (VNIR) hyperspectral imager to acquire multiple HSIs. Each HSI pixel included 462 wavelength bins ranging from 391 to 1011 nm. To evaluate the ability of GSM to identify potential contaminant sources, rhodamine dye, a commonly used tracer in hydrological studies, was released into the pond. Two UAV flights spaced 15 min apart were used to capture the evolution of the plume as it dispersed into the surrounding water.

The hyperspectral imager is in a pushbroom configuration so that each HSI was captured one scan line at a time. Data from an embedded GPS/INS unit enable direct georectification of captured imagery using the method described in [44]. Additionally, an upward-facing Ocean Optics UV-Vis NIR spectrometer with a cosine corrector was included on the top of the UAV to measure incident solar irradiance. The configuration of the UAV together with a sample HSI are shown in Figure 4. For a more detailed description of the system, the reader is directed to references [11,12].

Raw HSIs were converted to reflectance using the downwelling irradiance spectrum captured simultaneously with each HSI. Given the UAV flies with the imager oriented to nadir, the reflectance is then given by

$$\rho \left(\lambda \right)=\pi L\left(\lambda \right)/E_d\left(\lambda \right)$$

(20)

where L is the spectral radiance, $E_d$ is the downwelling irradiance, and a factor of $\pi$ is included resulting from the assumption of diffuse upwelling radiance [45]. The UAV flights were performed near solar noon to maximize the amount of sunlight illuminating the water. For this pond in North Texas, this corresponded to an average solar zenith angle of $56.7^{\circ }$, resulting in HSIs with negligible sunglint effects.

From the collected HSIs, a water-only pixel mask was generated by identifying all pixel spectra with a normalized difference water index (NDWI) greater than $0.25$ as defined in ref. [46]. Of these water pixels, a combined data set of $15,000$ spectra was sampled for GSM training. As a final processing step, reflectance spectra were limited to $\lambda \le 900$ nm, as wavelengths above this threshold showed significant noise.

In order to justify values for the appropriate number of endmembers, the final data set was decomposed using principal component analysis (PCA). A plot of the explained variance organized for each PCA component is shown in Figure 5.

The PCA decomposition of the data suggests that at least 3 endmembers should be used for a mixing model, and beyond 6 there is little added benefit. Based on these observations, multiple GSM models were trained with $N_v$ ranging from 3 to 6, ${\lambda }_e$ ranging from $0.001$ to $1.0$, and ${\lambda }_w$ ranging from 1 to 1000 in order to explore the GSM parameter space. From these, a final model was identified using the BIC, AIC, and reconstruction RMSE. The resulting GSM was then explored to examine extracted endmembers and map the evolution of the rhodamine plume by using the abundances given by the latent space representation of each pixel in the HSIs.

4. Results

4.1. Linear Mixing

The results of GSM training on the synthetic linear mixing data set described in Section 3 are illustrated in Figure 6. Three versions of NMF were trained corresponding to Euclidean (${\ell }_2$), KL-Divergence, and ${\ell }_{2,1}$ cost functions. GSM models using both a regular simplex grid and a grid of points sampled using a uniform Dirichlet distribution (referred to as a big GSM model) were trained to compare performance in linear mixing tasks.

The quality of endmember extraction is measured by the mean spectral angle and the mean endmember RMSE. As Figure 6 indicates, both versions of the GSM outperformed their NMF counterparts. Additionally, we note that for all GSM models, even including $\mathrm{SNR}=0$, all model weights $W_{dm}$ for $m>N_v$ corresponding to non-linear mixing were driven identically to $0.0$. This confirms that, for data sets with purely linear mixing and random noise, the GSM correctly fits a mixing model without introducing unnecessary complexity.

The quality of the unmixing, that is, the estimation of the abundance, performed by each model was evaluated using the mean abundance RMSE. Here, we again see that both versions of the GSM outperformed NMF. To justify that all models were fairly trained to convergence, the data reconstruction RMSE was also computed. This metric uses the trained mixing model to compute the error between the original data set and the reconstructed spectra generated using the extracted endmembers and their abundances. For all GSM and NMF models, the data reconstruction RMSE converged to the level of random noise introduced into the data. These results reflect a fair comparison between the NMF and GSM models.

In Figure 7, we plot the endmembers extracted for the GSM model trained on the synthetic data set with $\mathrm{SNR}=20$.

The extracted endmembers clearly fit the original source spectra while capturing local reflectance features. Furthermore, the GSM precision parameter $\beta$, which is tuned during model training, provides an assessment of spectral variability due to random noise. In Figure 7, this is indicated by the colored band centered around each extracted spectrum with a width of $2\sqrt{{\beta }^{-1}}$ corresponding to 2 standard deviations. The SNR of 20 added to this example corresponds to zero-mean Gaussian noise with a standard deviation of $\sigma =0.0493$. After training, the GSM found $\sqrt{{\beta }^{-1}}=0.0495$, accurately identifying this introduced noise. The ability to assess the spectral variability of extracted endmembers is a key advantage resulting from the probabilistic formulation of the GSM.

4.2. Non-Linear Mixing: Rhodamine Dye Plume

For the data set of real HSI spectra described in Section 3, 80 GSM models were trained to explore the GSM hyperparameter space. The performance of the model was compared using the BIC, AIC, and reconstruction RMSE with the results of the top 10 performing models shown in

Table 1. GSM hyperparameter optimization for the real HSI data set: Multiple GSM models were trained to explore the impact of model hyperparameters. $N_v$ values from 3 to 6 were explored as suggested by the PCA decomposition of the data set. ${\lambda }_e$ was varied from $0.001$ to $1.0$ and ${\lambda }_w$ ranged from 1 to 1000. Here, we report the top 10 models ranked by increasing BIC. The AIC and reconstruction RMSE are also included for comparison.

${\mathit{N}}_{\mathit{v}}$	${\mathit{\lambda }}_{\mathit{e}}$	${\mathit{\lambda }}_{\mathit{w}}$	BIC	AIC	Reconstruction RMSE
3	$0.01$	$1.0$	$-6.195\times {10}^7$	$-6.269\times {10}^7$	$0.000989$
3	$0.001$	$1.0$	$-6.194\times {10}^7$	$-6.268\times {10}^7$	$0.000989$
3	$0.1$	$1.0$	$-6.192\times {10}^7$	$-6.265\times {10}^7$	$0.000991$
3	$1.0$	$1.0$	$-6.186\times {10}^7$	$-6.260\times {10}^7$	$0.001190$
4	$1.0$	$1.0$	$-6.181\times {10}^7$	$-6.255\times {10}^7$	$0.001002$
4	$0.1$	$1.0$	$-6.175\times {10}^7$	$-6.249\times {10}^7$	$0.001008$
4	$0.01$	$1.0$	$-6.173\times {10}^7$	$-6.247\times {10}^7$	$0.001009$
4	$0.001$	$1.0$	$-6.173\times {10}^7$	$-6.247\times {10}^7$	$0.001009$
4	$0.1$	$10.0$	$-6.171\times {10}^7$	$-6.245\times {10}^7$	$0.001011$
3	$1.0$	$10.0$	$-6.166\times {10}^7$	$-6.239\times {10}^7$	$0.001014$

.

As the table indicates, a non-linear GSM with $N_v=3$, ${\lambda }_e=0.01$, and ${\lambda }_w=1.0$ achieved minimum BIC, AIC, and reconstruction RMSE values. The lower value of ${\lambda }_w$ identified from these models reflects the presence of non-linear mixing effects in the HSI data. Although the values of model weights $W_{dm}$ for $m>N_v$ corresponding to non-linear mixing were 0 for the synthetic data set, here these weights obtained a small but non-negligible median value of $0.0012$.

Reflectance spectra generated by the trained GSM corresponding to the maximum abundance values for each endmember are plotted in Figure 8. Based on their signatures, these endmembers are identified with water, vegetation (including near shore, filamentous blue-green algae), and the rhodamine dye plume.

By assigning a unique color to each endmember, the spatial distribution of HSI spectra can be visualized by using estimated abundances to smoothly interpolate between endmember colors as shown in panel b of Figure 8. In other words, this method provides a fuzzy semantic segmentation of HSI pixels. Alternatively, each HSI pixel could be assigned a single class corresponding to the endmember with maximal abundance. In this way, the GSM can be used to visualize high-dimensional HSI spectra.

To showcase the ability of the GSM to identify spatially localized contaminant sources, abundance maps were generated for each individual end member using the trained GTM to unmix all water pixels in the HSI. In Figure 9, these abundance maps are compared for each endmember. The abundance map for the water endmember covers a majority of the center of the pond and decreases near the edge of the rhodamine plume where the dye and water mix. The vegetation in shallow water near the shore and the rhodamine dye plume in the western half of the pond are also clearly identified by their abundance maps. Given an HSI pixel resolution of $0.1\times 0.1$${\mathrm{m}}^2$, the spatial extent of vegetation is estimated to be $378.6$ ${\mathrm{m}}^2$ while the plume initially extends across $255.7$ ${\mathrm{m}}^2$.

Applying the trained GSM to unmix HSI from UAV flights provides an efficient way to map the dispersion and transport of contaminant. Figure 10 demonstrates this by mapping the growing extent of rhodamine dye between successive UAV flights. In just 15 min, the plume expands from an initial area of $255.7$${\mathrm{m}}^2$ to over $571.8$${\mathrm{m}}^2$ as the dye mixes with the surrounding water.

5. Discussion

The proliferation of hyperspectral imaging technologies in the remote sensing and environmental monitoring communities underscores the need for efficient algorithms that can make sense of high-dimensional HSIs. Of particular interest are unsupervised algorithms, which utilize all available wavelength bins to identify unique source profiles that present realistic scenes. The GSM introduced in this paper is a novel method that simultaneously performs endmember extraction and spectral unmixing. Unlike popular endmember extraction algorithms, such as VCA, PPI, and N-FINDR, GSM does not assume the presence of pure pixels in HSIs. Furthermore, the flexible structure of the mapping from the GSM latent space to the spectral data space allows the GSM to model non-linear mixing effects, distinguishing it from other widely used linear models such as NMF. Being a probabilistic model also enables the GSM to quantify the spectral variability through the precision parameter $\beta$ as indicated in Figure 7. This additional capability is critical for analysis of realistic HSI where sensor noise, viewing geometry, and scene illumination all introduce uncertainty into extracted endmembers.

The growing popularity of deep learning methods has led many to explore applications of deep learning to the spectral unmixing problem. Variations in auto-encoder architectures are a popular choice for both endmember extraction and unmixing. For example, Borsoi et al. used a variational autoencoder to more accurately identify the spectral variability in an LMM [28]. Similarly, Palson et al. introduced an autoencoder with convolution layers to extend the LMM by incorporating spatial context from neighboring HSI pixels [29]. However, the increased complexity of deep neural network approaches makes interpreting the latent space learned by encoder layers challenging and can dramatically increase training time and data volume requirements. In contrast, the latent space of the GSM is immediately interpretable due to the imposed simplex geometry, while the form of $\psi$ leaves room to explore a variety of non-linear mixing models. In this paper, $\psi$ was constructed to add non-linear effects as additions or perturbations to an underlying linear mixing model, with the ${\lambda }_w$ hyperparameter provided to control the degree of non-linearity considered. In the original GTM formulation that inspired the GSM, the mapping $\psi$ strictly uses non-linear features generated by Gaussian RBFs. Alternatively, one could easily construct $\psi$ for specific non-linear models such as bilinear mixing by introducing features that depend on higher-order combinations of the latent space (abundance) coordinates. However, it is important to note that the EM algorithm formulated for GSM training is made possible due to the linear combination of activations $\varphi \left(z\right)$ with model weights $\mathbf{W}$. For mixing models that instead depend non-linearly on weights $\mathbf{W}$, an EM algorithm may require iterative non-linear solves during each M step.

The main limitation of the GSM is the curse of dimensionality encountered when generating a grid on the $(N_v-1)$-simplex for large numbers of endmembers, $N_v$. This can be mitigated by instead randomly sampling points within the interior of the simplex using a uniform Dirichlet distribution to obtain a predetermined number of nodes across the latent space. As the mixing coefficients ${\pi }_k$ are adapted during training, variability in spacing between nodes should not significantly affect the performance of the model. This was confirmed for the simulated data set as shown in Figure 6. In terms of computational efficiency, each expectation step involves $\mathcal{O}(K\times N)$ operations to update the entries of the responsibility matrix, leading to extended training times for considerably large data sets. This can be addressed by augmenting the EM procedure to use mini-batches of training samples, as outlined by Bishop et al. for the GTM in ref. [47]. Rather than updating the full responsibility matrix during each iteration, a subset of $\mathbf{R}$ corresponding to a single batch of training data can be evaluated with all other entries kept constant. The GSM may also be extended in other ways, for example, by replacing the precision parameter $\beta$ with a vector $\beta$ of parameters to model wavelength-dependent spectral variability common to many hyperspectral imaging platforms.

Based on the results of applying the GSM to real HSIs as illustrated in Figure 9 and Figure 10, one clear application of the GSM is contaminant identification and water quality assessment. Modern UAV platforms make it possible to rapidly image bodies of water where direct access for in situ data collection is restricted. By equipping UAVs with hyperspectral imagers, the GSM can be used to identify abnormal spectral signatures corresponding to localized contaminant sources. Generating semantic segmentations of collected HSIs as in panel b of Figure 8 makes it easy to compare multiple HSIs without resorting to pseudocolor images generated from a limited number of wavelength bands. Since the GSM models the distribution of all HSI spectra rather than individual parameters, it may also aid in situ data collection by suggesting sampling points that maximize the area traversed in the GSM latent space rather than uniformly sampling across a wide spatial extent. Similar approaches have been developed to guide autonomous data collection vehicles by casting route planning as a prize collecting traveling salesman problem subject to resource constraints [48,49].

As a final consideration, we note that the problem of endmember extraction and spectral unmixing for HSIs is identical to source apportionment in the context of air quality. Here, the goal is to identify measurement profiles associated with sources of ambient air pollution based on measurements at a receptor site. A popular model for source apportionment studies is Positive Matrix Factorization (PMF) introduced by Paatero and Tapper [50,51]. Just like NMF, PMF decomposes measurements into a convex combination of source profiles and relative abundances subject to non-negativity constraints. These linear receptor models assume that the sources are not transformed during transport to the receptor, ignoring changes due to chemical reactions. Non-linear mixing models such as the GSM introduced here may, therefore, prove beneficial in this additional domain.

6. Conclusions

In this paper, we introduced the GSM as a novel method for unsupervised endmember extraction and spectral unmixing. The model was compared to three varieties of NMF on a synthetic data set of linearly mixed spectra from the USGS spectral database. Across all noise levels, the GSM outperformed NMF for both endmember accuracy and abundance estimation. When applied to real HSIs collected over a north Texas pond, the GSM accurately models spectral mixing, including non-linear effects. Rhodamine dye was released into the pond to simulate a localized contaminant source. Endmember abundances estimated by the GSM corresponding to the dye accurately tracked the evolution of the plume as it mixed into the surrounding water. These examples illustrate the power of combining non-linear mixing models together with UAV-based hyperspectral imaging for environmental monitoring and water quality analysis. Future work will further develop the GSM to incorporate wavelength-dependent spectral variability. Additionally, the GSM may prove useful for air quality studies as a non-linear source apportionment model.

7. Patents

A provisional US patent application for this work has been filed on 2024-10-10 entitled Generative Simplex Mapping: Non-Linear Endmember Extraction And Spectral Unmixing For Hyperspectral Imagery, application number 63/705,854.