Spectral Variability Augmented Sparse Unmixing of Hyperspectral Images

The Linear Mixture Model (LMM), which is well-known for its simplicity and explicit physical meaning, has been widely utilized to model the relationship between mixtures and endmembers. In the LMM, the photons reflected from different ground objects within one pixel are assumed not to interfere with each other. As a result, the observed hyperspectral image matrix $\mathbf{R}\in\mathcal{R}^{b\times n}$ is assumed to be the product of endmember matrix $\mathbf{E}\in\mathcal{R}^{b\times p}$ and its corresponding abundance matrix $\mathbf{A}\in\mathcal{R}^{p\times n}$ plus a noise matrix $\mathbf{N}\in\mathcal{R}^{b\times n}$. Their relationship can be formulated as:

in which $b$ is the number of bands, $p$ is the number of endmembers, and $n$ is the number of pixels in the image. Here, the $i$-th column vectors of $\mathbf{E}$ and $\mathbf{A}$, denoted by $\mathbf{e}_{i}$ and $\mathbf{a}_{i}$, correspond to the $i$-th endmember and the abundance of the $i$-th spectral pixel $\mathbf{r}_{i}$, respectively.

In order to make the LMM physically meaningful, the abundance sum-to-one constraint (ASC) and the abundance nonnegative constraint (ANC) must be imposed, which can be expressed as:

Generally, the ASC can be embedded into this model by adding an additional pseudo band to the data matrix and the endmember matrix [53]. Therefore, in the following analysis, only the ANC is explicitly considered.

Due to the wide availability of spectral libraries, sparse unmixing models use mathematical relaxations that are computationally easier to solve. Let $\mathbf{M}\in\mathbb{R}^{b\times m}$ be a large spectral library containing $m$ spectral signatures, and $\mathbf{A}\in\mathbb{R}^{m\times n}$ denotes the abundance maps corresponding to library $\mathbf{M}$ for the observed data $\mathbf{R}$. As the number of endmembers involved in a mixed pixel is usually very small when compared with the size of the spectral library, the abundance matrix $\mathbf{A}$ contains many zero values, which means $\mathbf{A}$ is sparse. Thus, the unmixing problem can be formulated as an $l_{2}-l_{0}$ optimization problem,

where $\|\cdot\|_{F}$ is the Frobenius norm, $\|\cdot\|_{0}$ is the $l_{0}$-norm which denotes the number of nonzero components of a vector, and $\lambda$ is a regularization parameter used to adjust the weight of the sparsity in this model.

Problem (4) is nonconvex and difficult to be solved. However, the $l_{0}$-norm constrained sparsity can be relaxed using $l_{1}$ norm, by which the sparse unmixing algorithm via variable splitting and augmented Lagrangian (SUnSAL) [9] using the $l_{2}-l_{1}$ norm to solve the unmixing problem as follows,

where $\|\mathbf{A}\|_{1,1}=\sum_{j=1}^{n}\|\mathbf{a}_{j}\|$ with $\mathbf{a}_{j}$ being the $j$-th column of $\mathbf{A}$. However, the Problem (5) is heavily influenced by the high correlation of spectral libraries due to the underdetermined nature. To admit a sufficiently sparse solution for sparse unmixing will guarantee a more accurate abundance estimation. In [16], the Collaborative Sparse Unmixing via variable Splitting and Augmented Lagrangian (CLSUnSAL) algorithm imposes the joint sparsity with an $l_{2,1}$ mixed norm among the endmembers simultaneously for all of the pixels, whose objective function can be defined as follows,

where $\mathbf{a}^{i}$ denotes the $i$-th line of $\mathbf{A}$ and $\sum_{i=1}^{m}\|\mathbf{a}^{i}\|_{2}$ is the so-called $l_{2,1}$ mixed norm.

In a hyperspectral image, pixels from a homogeneous ground object often present in adjacent areas, which result in lots of pure spatial neighborhoods in the image. If all these pure spatial neighborhoods can be detected, their representative spectral signatures can be selected as endmember candidates. Different endmember candidates can be considered as coming from homogeneous ground objects provided that the spatial neighborhoods they represent cover pixels in common [54]. Based on this spatial refinement scheme, the remaining spatially independent endmember candidates become more representative since they represent larger areas. The number of endmember candidates to constrcut the in-situ spectral library will also be significantly reduced. Note that the in-situ spectral library $\mathbf{X}$ is extracted from the observed image $\mathbf{R}$ (refer to [54]), containing multiple spectra of each possible material (a total of $p$ materials), written as

in which $\mathbf{X}^{(i)}$ represents the library subset corresponding to the $i$-th possible material and is composed of $m$ spectral signatures of the current material.

Through the previous section, we have obtained the image spectral library including the spatially independent and spectrally pure pixels in the homogenous region according to the pixel purity. Principal Component Analysis (PCA) is a classic method widely used in signal feature extraction, data dimensionality reduction and compression. It retains the main information of the original image data by retaining the first $k$ principal components corresponding to the larger feature value, and the remaining principal components correspond to noise, spectral variability, and/or outliers.

In this section, PCA is performed on its covariance matrix after standardizing the pixel data $X$ of the in-situ spectral library as

where $\lambda_{1}>\lambda_{2}>\cdot\cdot\cdot>\lambda_{n}$ are the eigenvalues of its covariance matrix $\mathbf{C}$, and $\mathbf{W}$ is composed of eigenvectors of the corresponding eigenvalue matrix $\mathbf{D}$. Among all the principal components, the first $k$ correspond to the unique feature types of various pure pixels in the homogeneous area, while the other principal components correspond to the local spectral variability of the pixels in the homogeneous area.

The determination of $k$ can be achieved by comparing the size of the characteristic value. Set a threshold $\zeta$, if the cumulative eigenvalues corresponding to the first $k$ pivotal elements account for the percentage of the total eigenvalues to reach the threshold:

Therefore, the first $k$ principal components can be calculated by the first $k$ eigenvectors,

Then, the reconstructed spectral library by dominant characteristic components is termed as the endmember library,

Similarly, the next $n-k$ principal components are constructed as the spectral variability library as

Fig. 1 displays the flowchart of proposed SVASU where the well-known Jasper data set is taken as an example. The proposed SVASU adopts a two-order decomposition structure, in which the first order is to decompose the hyperspectral image data into the product of endmember library and abundance fractions roughly, the second order is to further represent the reconstruction error of pixels from first order as the influence of spectral variability library.

Let $\mathbf{V}\in\mathbb{R}^{b\times l}$ be the spectral variability library containing $l$ spectral variability signatures, and $\mathbf{B}\in\mathbb{R}^{l\times n}$ denotes the coefficients corresponding to library $\mathbf{V}$. Based on the objective function in Eq. (6), an extra spectral variability augmented data fitting term is introduced in the proposed SVASU model, along with the coefficient generalization regularized term, which is expressed as follows,

where $\alpha$ balances the contribution of data fitting from both spectral spectra and spectral variability, $\beta$ and $\gamma$ controls the sparsity of abundance and the generalization of spectral variability coefficient.

Following [55], we relax the term $\|\mathbf{A}\|_{2,1}$ by $\rm{Tr}(\mathbf{A}^{T}\mathbf{D}\mathbf{A})$. Thus the function in (13) can be written as

where $\mathbf{D}\in\mathbb{R}^{m\times m}$ is a diagonal matrix with the $i$-th diagonal element formulated as follows, in which $\epsilon>0$ is a stabilization parameter,

In order to integrate non-negative constraint conditions into the objective function $\bf{\Theta}$, we introduce Lagrangian multipliers $\delta\in\mathbb{R}_{+}^{m\times n}$ to restrict $\mathbf{A}\geq 0$.Therefore, the function (14) is equivalent to the following function:

When $\mathbf{D}$ is fixed, the original problem can be turned into a general convex optimization problem. To update all matrixes of the objective function, we use an alternating projected gradient method. By setting the partial derivatives of with respect to $\mathbf{A}$ and $\mathbf{B}$ respectively, we obtain the following function:

We take advantage of KKT conditions in this mathematical optimization problem. Therefore, the following condition should be satisfied:

Then, we design the following update rules. The convergence of two update rules will be analyzed in the APPENDIX.

Therefore, an iterative procedure is adopted. In each iteration, $\mathbf{A}$ is calculated with the current $\mathbf{D}$ and $\mathbf{B}$, and then $\mathbf{B}$ and $\mathbf{D}$ are updated based on the current calculated $\mathbf{A}$. The iterative procedure is repeated until the algorithm converges. The pseudo code of SVASU is illustrated in Algorithm 1.

To quantitatively evaluate the performance of our proposed SVASU algorithm, the two types of pixel reconstruction error are employed as in the experiments. Firstly, the root mean square error (RMSE) is used to measure the distance between the true signal and its estimated value, like that of image pixel $\rm{\mathbf{r}}$ and its reconstructed version $\rm{\hat{\mathbf{r}}}$ as $\rm{RMSE_{R}}$, and true abundance $\rm{\mathbf{a}}$ and its estimated version $\rm{\hat{\mathbf{a}}}$ as $\rm{RMSE_{A}}$, which are defined as,

Furthermore, the signal-to-reconstruction-error (SRE) is also adopted to measure the quality of the reconstruction of signals by different algorithms, which can demonstrate more information regarding the power of the error in relation to the power of the signal. Generally, the higher the SRE (dB) value, the better the performance of the algorithm. Similarly, we definite $\rm{SRE_{R}}$ and $\rm{SRE_{A}}$ as follows,

In this section, the simulated data are generated in the way similar to [56] and [57] based on the LMM in (1). The synthetic data set has $n=10000$ pixels. The endmember spectral signatures are using five spectra of minerals from the U.S. Geological Survey (USGS) containing as many as 420 bands covering from 400 to 2500 nm. Then, 30 dB Gaussian noise is added to each signature to simulate spectral variability. Considering that each pixel is unlikely to have a large number of endmembers in real scenarios (typically, less than 5), the maximum number of active endmembers per pixel is set to 4 to ensure abundance sparsity. Furthermore, the abundance vectors corresponding to active endmembers are generated using a mixture of uniform Dirichlet distributions [58]. Finally, 40 dB Gaussian noise is added to the whole image. In this experiment, the In-situ spectral library is simply composed of the selected five spectra from USGS and their corresponding spectra owning spectral variability. After the library segmentation by PCA, we acquire the (a) endmember library and (b) spectral variability library as shown in Fig. 2.

First of all, the impact of parameters, including $\alpha$, $\beta$, and $\gamma$, are discussed through several sets of simulated experiments. The curves of these parameters of our method are shown in Fig. 3(a) and (b). The evaluation metrics adopted here are the average SRE and RMSE calculated by (20)-(23). On one hand, in Fig. 3(a), $\alpha$ is preset to 10 and $\beta$ and $\gamma$ vary in $[1e0,1e3]$ and $[1e2,1e5]$ exponentially. The proposed SVASU algorithm achieves a desired result when $\beta=10$ and $\gamma=1e4$, with lower RMSE value and higher SRE value both for abundance and pixels. On the other hand, Fig. 3(b) displays the performances of SRE and RMSE under different values of parameter $\alpha$ tend to prominent when $\alpha$ reaches 9. Thus, the parameter $\alpha$, $\beta$, and $\gamma$ are respectively set to be 9, 10, and $1e4$ in the experiments on synthetic data set.

Table I shows the average SRE and RMSE value between the ground-truth pixels or abundance and corresponding matrixes acquired by different algorithms over synthetic data set. The performance of the proposed SVASU algorithm outperforms other methods on the estimation of either fractional abundance or pixel reconstruction.

The real-world Jasper data set is collected by an airborne visible/infrared imaging spectrometer (AVIRIS) sensor¹¹1Data available online at https://rslab.ut.ac.ir/data., which consists of 224 spectral bands ranging from 380 to 2500 nm and the spectral resolution is up to 9.46 nm. The size of this scene is 512$\times$614. In the experiments, a subimage of 100$\times$100 pixels is used. Due to dense water vapor and atmospheric effects, there exists some noisy bands, including the bands 1$-$3, 108$-$112, 154$-$166, and 220$-$224. After removing these bands, a total of 198 reflectance bands are finally adopted. There are four endmembers in the selected region, including road, soil, water, and tree, and their corresponding ground-truth abundance are shown in Fig. 4.

As shown in Fig. 1, the In-situ spectral library obtained by SPEE is shown in the first column, and then after the segmentation step described in Section III. B, we get the endmember library and spectral variability library shown in Fig. 5. To analyze the impact of three regularization parameters, including $\alpha$, $\beta$, and $\gamma$, we first fix $\alpha$ at 10 to discuss the performance of SVASU under different values of $\beta$ and $\gamma$. As shown in Fig. X(a), with the increase in $\beta$ value, the performances of abundance estimation and pixel reconstruction demonstrate different trends of development. Since the gap between the maximum and minimum values of evaluation matrixes for pixel reconstruction is small, the optimal parameter $\beta$ is set to $1e3$ and $\gamma$ is set to $1e6$ corresponds to good capability of abundance estimation. After that, we further discuss the range of $\alpha$ as shown in Fig. X(b). Through observation, it can be found that when the value of $\alpha$ is lower than 3, the unmixing results of both abundance accuracy and data fidelity are superior to other cases. Thus, the optimal value of $\alpha$ is set to 1.

Table II shows the average SRE and RMSE values between the ground-truth pixels or abundance and corresponding matrixes acquired by different algorithms over synthetic data set. The performance of the proposed SVASU algorithm outperforms other methods on the estimation of fractional abundance. However, due to the second level decomposition by our extracted spectral variability, the performance of pixel reconstruction is not as well as other sparse unmixing methods as expected.

Furthermore, to carefully analyze the results of two-order reconstruction, Fig. 6 depicts the visual results of the important items. Each visual result is averaged over each band, including the image data $\mathbf{R}$, the first order reconstruction item $\mathbf{MA}$, the second order reconstruction item $\mathbf{VB}$, the first order reconstruction error item $\mathbf{R-MA}$, and the second order reconstruction error item $\mathbf{R-MA-VB}$. It can be observed that on the Jasper data set, the amount of information occupied by $\mathbf{MA}$ is probably twice as much information as that represents by $\mathbf{VB}$. Moreover, through the second order reconstruction by spectral variability library, the image reconstruction error reduced to a tenth on Jasper data set, which completely proves the validity and rationality of the proposed SVASU method.

The hyperspectral data set named Cuprite data set, acquired by AVIRIS²²2http://aviris.jpl.nasa.gov/html/aviris.freedata.html., is also adopted in real experiments, which contains 224 bands ranging from 370 to 2510 nm with a ground instantaneous field of view of 20 m. The cropped image corresponds to a $350\times 350$ pixel subset of the sector labeled as $f970619t01p02_{r}02_{s}c03.a.rfl$ in the online data. After removing noisy bands and water-absorption bands (including bands $1-4$, $105-115$, $150-170$, and $223-224$), a total of 186 reflectance bands are finally adopted. The optimal parameters of proposed SVASU method are $\alpha=3$, $\beta=1e4$, and $\gamma=1e6$.

Given that there is no available ground-truth in terms of abundances and endmembers of Cuprite data set, only visual results of two-order reconstruction in SVASU are shown as Fig. 7. Similarly, each visual result is averaged over each band, including the two-order reconstruction item $\mathbf{MA}$, $\mathbf{VB}$, and the two-order reconstruction error item $\mathbf{R-MA}$, $\mathbf{R-MA-VB}$. It can be observed that the amount of information occupied by $\mathbf{MA}$ is more than twice as much information as that by $\mathbf{VB}$. This also confirms that in the processing of $PCA$, the main components are constructed by endmember library while the other components are constructed as the spectral variability library. Moreover, through the second order reconstruction by spectral variability library, the image reconstruction error reduced to less than a tenth on Cuprite data set. Therefore, the effectiveness of the proposed SVASU method are confirmed on the real-world data set.