Joint Sparse Sub-Pixel Mapping Model with Endmember Variability for Remotely Sensed Imagery

Abstract

Spectral unmixing and sub-pixel mapping have been used to estimate the proportion and spatial distribution of the different land-cover classes in mixed pixels at a sub-pixel scale. In the past decades, several algorithms were proposed in both categories; however, these two techniques are generally regarded as independent procedures, with most sub-pixel mapping methods using abundance maps generated by spectral unmixing techniques. It should be noted that the utilized abundance map has a strong impact on the performance of the subsequent sub-pixel mapping process. Recently, we built a novel sub-pixel mapping model in combination with the linear spectral mixture model. Therefore, a joint sub-pixel mapping model was established that connects an original (coarser resolution) remotely sensed image with the final sub-pixel result directly. However, this approach focuses on incorporating the spectral information contained in the original image without addressing the spectral endmember variability resulting from variable illumination and environmental conditions. To address this important issue, in this paper we designed a new joint sparse sub-pixel mapping method under the assumption that various representative spectra for each endmember are known a priori and available in a library. In addition, the total variation (TV) regularization was also adopted to exploit the spatial information. The proposed approach was experimentally evaluated using both synthetic and real hyperspectral images, and the obtained results demonstrate that the method can achieve better results by considering the impact of endmember variability when compared with other sub-pixel mapping methods.

1. Introduction

Mixed pixels are frequent in remotely sensed images due to coarse spatial resolution. Given the assumption that two or more different land-cover classes may exist in a single pixel, spectral unmixing techniques have been commonly used to quantitatively estimate the proportion of each pure spectral constituent (endmember) in a mixed pixel by establishing the relationship between the measured spectra of the pixel and the corresponding representative spectra of the land-cover classes [1,2]. The obtained abundance map can indicate the quantitative proportions of the land-cover classes in the mixed pixel, whereas the spatial distribution of each endmember in the pixel can be further estimated using sub-pixel mapping techniques.

Sub-pixel mapping, which in most cases can be regarded as a subsequent procedure to spectral unmixing, aims to arrange the spatial location of possible classes inside each mixed pixel given the obtained abundance map [3,4]. Under the principle of spatial dependence inspired by Tobler’s first law [5], sub-pixel mapping techniques can convert a low-resolution (LR) abundance map into a finer classification result by dividing the pixel into numerous smaller sub-pixels and allocating each sub-pixel to a certain class.

In the literature, many efforts have been directed towards the development of sub-pixel mapping techniques aimed at obtaining a finer classification map from a lower spatial resolution image [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23], such as the pixel-swapping algorithm [6,7], the Hopfield neural network [8,9], the spatial attraction model [10,11,12], genetic algorithms [13,14], multi-agent systems [15], maximum a posteriori (MAP) based techniques [16,17,18], and differential evolution [19]. However, most of these algorithms attempt to retrieve the finer map from a previously estimated set of abundance maps, which is commonly obtained by using spectral unmixing techniques. Therefore, the final sub-pixel mapping results strongly rely on the unmixing step, which is directly responsible for the quality of the abundance maps. While the two steps are performed independently, errors in the first step may propagate to the second one during the process. Moreover, an assumption is implicitly embedded in the “unmixing-mapping” procedure that the abundance map with higher accuracy always generates a better sub-pixel mapping result. However, this deduction needs to be further verified and therefore it is meaningful to develop sub-pixel mapping algorithms which act on the original remotely sensed images directly.

To address the inaccuracy of the utilized abundance map and improve the sub-pixel mapping performance, many researchers have attempted to incorporate the information of the original coarse remotely sensed image into the sub-pixel mapping model [24,25,26,27,28,29,30,31,32,33]. A Markov random field (MRF) model based approach was firstly introduced to generate super-resolution land cover maps from remote sensing data by Kasetkasem [24]. Then a supervised fuzzy c-means approach was proposed in [26] which incorporates the pixel-unmixing criterion in the objective function. To ensure the spectral consistency between the final sub-pixel mapping result and the original utilized coarse remotely sensed image, a spectral term based on the Gaussian maximum likelihood was incorporated into the sub-pixel mapping model [28]. In [29] the object boundaries, which are obtained from remotely sensed images by segmentation, are integrated into the sub-pixel model as an additional structural constraint. Additionally, a spectral difference between the original and the inversed signatures is calculated as part of the objective function of the commonly used genetic algorithm [30].

Although a better sub-pixel mapping result can be obtained by the aforementioned methods, which consider the inaccuracy of the abundance map, the abundance map is still needed as the intermediary part of the process for most methods. Therefore, it is essential to develop sub-pixel mapping algorithms which act on the original remotely sensed images directly. However, when it comes to the way of generating the sub-pixel mapping result from the original remotely sensed images directly, another issue arises which is a common problem in remote sensing researches. Owing to the environmental, atmospheric, and temporal factors, the spectral signatures of an identical material can vary greatly [34,35]. By ignoring these variations, errors are introduced and propagated throughout image analysis such as spectral unmixing and sub-pixel mapping. Therefore, it is necessary to consider the impact of endmember variability in remotely sensed image processing and many different methods that account for spectral variability have been studied such as the multiple endmember spectral mixture analysis [36], endmember bundles approaches [37], sparse unmixing [38] and multivariate statistical distribution based methods [39].

In this paper, a joint sparse sub-pixel mapping method (JSSM) was designed to combine the procedures of spectral unmixing and sub-pixel mapping together by introducing a so-called sub-pixel abundance map, which indicates the proportions of sub-pixels belonging to different land cover classes. In this way, the original remotely sensed image and the final sub-pixel map can be connected without the need for intermediate abundance maps, and the propagation of errors in the model can be mitigated. Moreover, the variable representative spectra for each endmember are selected as a priori and available in a library. In this way, endmember variability can be incorporated as a sparse constraint into the sub-pixel mapping model. Additionally, the total variation (TV) regularizer was also adopted to exploit the spatial characteristic of the deserved result. Since this model is built, the sub-pixel abundance map can be updated by the alternating direction method of multipliers (ADMM) algorithm and furthermore, the sub-pixel mapping result can be obtained.

The remainder of the paper is organized as follows. Section 2 outlines the sub-pixel mapping problem and Section 3 provides a detailed description of the proposed JSSM method. Section 4 explores the performance of the method using both synthetic and real hyperspectral data. Section 5 concludes the paper with some remarks and hints at plausible future research lines.

2. Sub-Pixel Mapping Problem

The key issue in a sub-pixel mapping problem is how to determine an optimal sub-pixel distribution of each class within a pixel. Spatial dependence has been regarded as a criterion that is related to the tendency of spatially close observations to be more alike than more distant observations [3]. As illustrated in Figure 1, given an abundance map obtained by spectral unmixing techniques, each coarse pixel can be divided into s × s sub-pixels, where s represents the scale factor. The number of sub-pixels for each land-cover class can be then determined by the fractional values of different classes. Figure 1 shows a sub-pixel mapping example with three classes. As shown in Figure 1a, a coarse pixel is divided into 16 (4 × 4) sub-pixels, where the scale fraction s equals 4, and 0.5 in the fraction image in red, which means that 8 (16 × 0.5) sub-pixels belong to land-cover class 1. Figure 1b,c describe two possible distributions of sub-pixels. Given the principle of spatial dependence, the former is perceived to be more optimal.

3. The Proposed Joint Sparse Subpixel Mapping Model

3.1. Basic JSSM Model

First of all, we define the formulations and notations used in this paper. Let $Y=\left[y_1,\cdot \cdot \cdot ,y_n\right]\in {\Re }^{b\times n}$ be the observed remotely sensed image with b spectral bands and n pixels, let $M=\left[m_1,\cdot \cdot \cdot ,m_p\right]\in {\Re }^{b\times p}$ collect the spectral signatures of p endmembers (p equals to the number of land cover classes if the endmember variability is not considered) present in $Y$, let $A=\left[a_1,\cdot \cdot \cdot ,a_n\right]\in {\Re }^{p\times n}$ stand for the abundance maps associating with $Y$ and $M$. Following the linear mixture model, we have

$$\begin{array}{l}Y=MA+N\\ s.t.:A\ge 0,1_p^{T}A=1_n^T\end{array}$$

(1)

where $N\in {\Re }^{b\times n}$ is the noise in the data, $A\ge 0$ and $1_p^{T}A=1_n^T$ are the so called nonnegative and sum-to-one constraints, $1_p^T={\left[1,1,\cdot \cdot \cdot ,1\right]}^T$ is a column vector of size p of 1s. Commonly, the endmember spectra are collected from the observed remotely sensed image manually or determined by endmember extraction methods such as vertex component analysis (VCA) [40] and N-FINDR [41].

Let $Z=\left[Z_1,\cdot \cdot \cdot ,Z_v\right]\in {\Re }^{p\times v}$ be the sub-pixel abundance map, where $v=n\times s^2$ is the number of sub-pixels in, $Z$ and s is the scaling factor. An explicit relationship between the abundance maps and the introduced sub-pixel abundance map can be established by resorting to a downsampling matrix as follows:

$$A=ZW,$$

(2)

where $W=\left[d_1,\cdot \cdot \cdot ,d_n\right]\in {\Re }^{v\times n}$ is the downsampling matrix, which can be constructed as

$$W=\frac{{(I_l\otimes 1_s^T)}^T\otimes {(I_r\otimes 1_s^T)}^T}{s^2},$$

(3)

where l and r are the number of lines and columns in $Y$, respectively, and the total number of samples in the observed image $Y$ is $n=l\times r$. In (3), $\otimes$ denotes the Kronecker operator, $I$ is an identity matrix with suitable dimension.

Then a joint model can be generated by combining (1) and (2):

$$\begin{array}{l}Y=MZW+N\\ s.t.:Z\ge 0,1_p^{T}Z=1_v^T\end{array}$$

(4)

where $Z\ge 0$ and $1_p^{T}Z=1_v^T$, similar to those in (1), are the nonnegative and sum-to-one constraints, respectively. In this way, we can directly associate the final sub-pixel abundance map $Z$ with the coarser image $Y$, along with the constraints.

As mentioned before, it is essential to take the endmember variability into consideration for the joint sub-pixel mapping model. In this paper, the popular sparse representation theory was utilized. The idea of sparse representation was firstly introduced for spectral unmixing of hyperspectral images considering that pure signatures may not present in the input data [42,43]. It assumes that mixed pixels can be expressed in the form of linear combinations of a number of pure spectral signatures known in advance and available in a library. In this proposed JSSM model, a number of endmember spectra are selected from the remotely sensed images manually for each land cover class to represent the endmember variability and an endmember library can therefore be generated. Then the fixed endmember matrix $M=\left[m_1,\cdot \cdot \cdot ,m_p\right]\in {\Re }^{b\times p}$ in (4) is replaced by the spectral library $D=\left[d_1,\cdot \cdot \cdot d_q\right]\in {\Re }^{b\times q}$ as (5) shows. Specifically, $q={\displaystyle \sum _{i=1}^p{q}_i}$ is the number of endmember spectra in the library in which p is the observed number of land cover classes and $q_i$ is the number of collected spectra for certain class i.

$$Y=DZW+N$$

(5)

Due to the small number of endmembers contributing to a mixed pixel, $Z\in {\Re }^{q\times v}$ is extremely sparse due to most of the lines being full of zeros. This sparse problem for the hyperspectral image can therefore be written as:

$$\underset{Z}{\min }{\Vert Z\Vert }_{0}s.t.{\Vert Y-DZW\Vert }_F^2\le \delta$$

(6)

where ${\Vert Z\Vert }_0$ denotes the number of nonzero components of $Z$ and $\delta \ge 0$ is the tolerated error derived from the noise or the model itself. Commonly the $L_0$ term can be replaced by the $L_1$ norm and this problem can be restated as follows [44,45]:

$$\underset{Z}{\min }{\Vert Z\Vert }_{1}s.t.{\Vert Y-DZW\Vert }_F^2\le \delta$$

(7)

where ${\Vert Z\Vert }_1={\displaystyle \sum \left|Z_{i,j}\right|}$ and the above constraint optimization problem can be converted into an unconstrained version by minimizing the respective Lagrangian function:

$$\underset{Z}{\min }\frac{1}{2}{\Vert Y-DZW\Vert }_F^2+\lambda {\Vert Z\Vert }_1+{\tau }_{R+}(Z)$$

(8)

where the parameter $\lambda$ is a non-negative parameter which controls the relative weight of the sparsity of the solution and the term ${\tau }_{R+}(Z)$ is the non-negative constraint. Moreover, the sum-to-one constraint is not considered as suggested by [42,43,46,47]. The term ${\tau }_{R+}(Z)$ is defined in the set $R+$ and obeys rules such as ${\tau }_{R+}(Z)=0$ when $Z\in R+$, otherwise ${\tau }_{R+}(Z)=+\infty$.

3.2. Spatial Prior Constraint

Owing to the ill-posed nature of the sub-pixel mapping problem which aims to generate a fine map from the coarser image, a spatial prior constraint is further imposed to regularize the problem. In this paper, the anisotropic TV model [48] is utilized by promoting piecewise smooth transitions in the sub-pixel abundance map. Moreover, it would be more suitable to conduct this operator on the class domain and not the endmember domain considering its physical interpretation as Figure 2 shows. It can be seen that the sub-pixel abundance map for all endmembers in the library should be first grouped according to the categories of endmembers for different classes.

Then the utilized TV operator can be imposed on Z* as (9).

$$TV(\mathrm{Z*})={\Vert {\nabla }_x\mathrm{Z*}\Vert }_1+{\Vert {\nabla }_y\mathrm{Z*}\Vert }_1,$$

(9)

where ${\nabla }_x$ and ${\nabla }_y$ are linear operators denoting the horizontal and vertical first-order differences. It is obvious that only the two directions are taken into consideration to represent the spatial dependence. For pixel i of class k in image Z*, its two-dimensional position can be indicated as (m, u) in which $i=m\times n\times s+u$ and $i\le v$. Then each elements in ${\nabla }_x\mathrm{Z*}$ and ${\nabla }_y\mathrm{Z*}$ can be computed as $\nabla Z_{k,i}^x=Z_k[m+1,u]-Z_k[m,u]$ and $\nabla Z_{k,i}^y=Z_k[m,u+1]-Z_k[m,u]$, respectively, where $m=floor(i/(n\times s))$ and $u=i-m\times n\times s$. By combing the two linear operators, this TV model can be defined as:

$$TV(\mathrm{Z*})={\Vert \begin{array}{l}{\nabla }_x\mathrm{Z*}\\ {\nabla }_y\mathrm{Z*}\end{array}\Vert }_1={\Vert \nabla \mathrm{Z*}\Vert }_1={\Vert \nabla GZ\Vert }_1$$

(10)

3.3. Optimization

With these definitions, this optimization problem can be rewritten as

$$\underset{Z}{\min }\frac{1}{2}{\Vert Y-DZW\Vert }_F^2+{\lambda }_{TV}{\Vert \nabla GZ\Vert }_1+\lambda {\Vert Z\Vert }_1+{\tau }_{R+}(Z)$$

(11)

where ${\lambda }_{TV}$ is the regularization parameter for the TV term. To solve this problem, the widely used alternating direction method of multipliers [43,49] is employed to decompose (11) into a sequence of simpler ones. By introducing auxiliary variables, the following constrained problem can be obtained.

$$\begin{array}{l}\underset{Z,V_1,V_2,V_3,V_4,V_5,V_6}{\min }\frac{1}{2}{\Vert Y-D{V}_2\Vert }_F^2+{\lambda }_{TV}{\Vert V_4\Vert }_1+\lambda {\Vert V_5\Vert }_1+{\tau }_{R+}(V_6),\\ \begin{array}{ll}subjectto& V_1=Z\\ & V_2=V_{1}W\\ & V_3=GZ\\ & V_4=\nabla V_3\\ & V_5=Z\\ & V_6=Z\end{array}\end{array}$$

(12)

The expanded augmented Lagrangian formula is then obtained in (13) by adding the many ancillary terms that utilize the ADMM in the optimization problem (12).

$$\begin{array}{l}J(Z,V_1,V_2,V_3,V_4,V_5,V_6,S_1,S_2,S_3,S_4,S_5,S_6)\\ =\frac{1}{2}{\Vert Y-D{V}_2\Vert }_F^2+{\lambda }_{TV}{\Vert V_4\Vert }_{1,1}+\lambda {\Vert V_5\Vert }_{1,1}+{\tau }_{R+}(V_6)\\ +\frac{\mu }{2}\left(\begin{array}{l}{\Vert Z-V_1-S_1\Vert }_F^2+{\Vert V_{1}W-V_2-S_2\Vert }_F^2+{\Vert GZ-V_3-S_3\Vert }_F^2\\ +{\Vert \nabla V_3-V_4-S_4\Vert }_F^2+{\Vert Z-V_5-S_5\Vert }_F^2+{\Vert Z-V_6-S_6\Vert }_F^2\end{array}\right)\end{array}$$

(13)

where $\mu$ is the penalty parameter. In the following part, we elaborate on how to split the complicated optimization relation (13) based on the ADMM strategy, and give details of the derivation of the JSSM.

To achieve the solution Z of the augmented Lagrangian formula by performing an alternating minimization process, the partial derivative of (13) can be written as:

$$\begin{array}{l}{Z}^{k+1}\leftarrow \underset{Z}{\min }J(Z,V_1^k,V_2^k,V_3^k,V_4^k,V_5^k,V_6^k,S_1^k,S_2^k,S_3^k,S_4^k,S_5^k,S_6^k)\\ \Rightarrow Z^{k+1}\leftarrow {(G^{T}G+3I)}^{-1}*\left\{(V_1^k+S_1^k)+G^T(V_3^k+S_3^k)+(V_5^k+S_5^k)+(V_6^k+S_6^k)\right\}\end{array}$$

(14)

The optimization of $V_i$ can be computed as follows, which is obtained by means of derivation. First, we compute the optimization of $V_1^{k+1}$ as follows:

$$V_1^{k+1}\leftarrow \left\{(Z^{k+1}-S_1^k)+(V_2^k+S_2^k)W^T\right\}{(I+W{W}^T)}^{-1}$$

(15)

Similarly, the optimization of $V_2^{k+1}$ and $V_3^{k+1}$ can be written as (16) and (17).

$$V_2^{k+1}\leftarrow {(D^{T}D+\mu I)}^{-1}(D^{T}Y+\mu (V_1^{k+1}W-S_2^k))$$

(16)

$$V_3^{k+1}\leftarrow {({\nabla }^T\nabla +I)}^{-1}(G{Z}^{k+1}-S_3^k+{\nabla }^T(V_4^k+S_4^k))$$

(17)

Here, $\nabla$ is a convolution operator and it can be computed using the discrete Fourier transform diagonalization band by band as in [43]. Moreover, as the $L_1$ norm is not differentiable, in the case of $L_1$ sparse approximation of $V_4$ and $V_5$, the solution is obtained by a soft shrinkage operator, as follows:

$$V_4^{k+1}\leftarrow soft(\nabla V_3^{k+1}-S_4^k,\frac{{\lambda }_{TV}}{\mu })$$

(18)

$$V_5^{k+1}\leftarrow soft(Z^{k+1}-S_5^k,\frac{\lambda }{\mu })$$

(19)

According to the definition of ${\tau }_{R+}$, the optimization of $V_6$ can be:

$$V_6^{k+1}\leftarrow \max (Z^{k+1}-S_6^k,0)$$

(20)

Finally, the Lagrangian multipliers are updated:

$$\begin{array}{l}{S}_1^{k+1}\leftarrow S_1^k-Z^{k+1}+V_1^{k+1}\\ S_2^{k+1}\leftarrow S_2^k-V_1^{k+1}W+V_2^{k+1}\\ S_3^{k+1}\leftarrow S_3^k-G{Z}^{k+1}+V_3^{k+1}\\ S_4^{k+1}\leftarrow S_4^k-\nabla V_3^{k+1}+V_4^{k+1}\\ S_5^{k+1}\leftarrow S_5^k-Z^{k+1}+V_5^{k+1}\\ S_6^{k+1}\leftarrow S_6^k-Z^{k+1}+V_6^{k+1}\end{array}$$

(21)

In summary, the ADMM strategy can be utilized as indicated in Equations (14)–(21) iteratively until a given number of iterations or a tolerance value of adjacent results is achieved. Moreover, the ADMM method is initialized as $Z^0=V_i^0=S_i^0=0,i\in [1,6]$ and the maximum iterative time is set as 200 in case the given stop threshold cannot be satisfied.

As a result, the sub-pixel abundance map Z can be obtained and it should be further grouped as Figure 1 shows to obtain the sub-pixel abundance map Z* in class domain where the value $Z_{k,i}^{*}$ in Z* denotes the proportion of pixel i belonging to class k. Then the sub-pixel mapping result $O$ with p classes can be expressed as $O=f(\mathrm{Z*})$, which means that it can be generated with the sub-pixel abundance map Z*, and different functions $f(\cdot )$ can be used to convert the abundance map to a final sub-pixel mapping result. In this paper, a simple winner-takes-all strategy was designed to generate the final sub-pixel mapping result $O$ as follows,

$$O_i=jif{Z}_{j,i}^{*}=\max \left\{Z_{k,i}^{*}|k\in (1,p)\right\}$$

(22)

4. Experiments and Analysis

The proposed JSSM has been compared with three traditional sub-pixel mapping algorithms, including the attraction model (AM) [10], Hopfield neural network (HNN) [8], and MRF [24], in which AM and HNN are based on an obtained abundance map while the MRF can be conducted on the original remotely sensed image directly. Additionally, the abundance fractions in the abundance map act as a mandatory constraint for the sub-pixel mapping result generated by the AM method, while for HNN, a soft constraint strategy is utilized which means the final sub-pixel mapping result does not need to satisfy the abundance fraction constraint. Moreover, the abundance map can also be classified and further rescaled to generate a finer classification map which was also used to compare with the proposed JSSM method and this strategy is named as the hard classification (HC) method. Additionally, two spectral unmixing methods, FCLS [50] and SUnSAL-TV [43], were employed to generate the abundance maps for AM, HC, and HNN. Moreover, the endmember library used in the JSSM method was also applied for SUnSAL-TV while for FCLS, only one representative spectrum was determined for each endmember. Owing to the fact that the spectral curves of the identical land cover class vary greatly and it is hard to determine a most representative endmember for each class, the average spectral signature for each class is generate as the input of FCLS. For the utilized HNN method, the parameters were set as follows: λ is 100, k1 = k2 = k3 = 1, dt = 0.01. Moreover, different parameters were tested for SUnSAL-TV and the abundance map, with the highest accuracy on comparing with the reference abundance map, was selected as the one used for the following sub-pixel mapping procedure. In the proposed JSSM method, three parameters were taken into consideration as mentioned above: the spatial constraint parameter ${\lambda }_{TV}$, the sparsity constraint parameter $\lambda$, and the penalty parameter $\mu$.

In our experiments, three synthetic and one real hyperspectral images were used to evaluate the proposed JSSM method in comparison to other techniques. For the synthetic image experiments, an original high-resolution remotely sensed image was firstly degraded to obtain a low-resolution image by applying an averaging filter. The low-resolution image was then used to obtain different sub-pixel mapping results with spectral unmixing and sub-pixel mapping techniques. The high-resolution remotely sensed image can be classified to generate a reference classification map to evaluate different sub-pixel mapping methods. For the real experiment, spectral unmixing and sub-pixel mapping algorithms were applied on the low resolution hyperspectral image and the classification result of a high resolution image, which covers an identical area to the low resolution one, was used as the reference map.

Our accuracy assessment was undertaken using the overall (OA), average (AA) and individual classification (IA) accuracies, as well as the Kappa coefficient. Moreover, the root-mean-square error (RMSE) index was also utilized to evaluate the generated abundance results with two different spectral unmixing techniques by comparing with the reference abundance map which was generated from the reference classification map.

In addition, to test the statistical significance of the differences in accuracy for results of the proposed method and other algorithms, the McNemar’s test [51] was used to compare the misclassification rates with the different methods. For the two classification maps $C_1$ and $C_2$, the McNemar’s test compares the number of pixels misclassified in $C_1$, but not in $C_2$($M_{12}$), with the number of pixels misclassified in $C_2$ while not in $C_1$($M_{21}$). If $M_{12}+M_{21}\ge 20$, the $X^2$ can be considered as a chi-squared distribution (with one degree of freedom) [52], [53] as (23):

$$X^2=\frac{{(\left|M_{12}-M_{21}\right|-1)}^2}{M_{12}+M_{21}}\approx {\chi }_1^2$$

(23)

The McNemar’s test accepts the hypothesis that the two classification methods have the same error rate at significance level $\epsilon$ if the value is less than or equal to ${\chi }_{\epsilon ,1}^2$ [54]. In other words, if the McNemar’s value is greater than ${\chi }_{\epsilon ,1}^2$, the two classification algorithms are significantly different. In this paper, the significance level $\epsilon$ is set as 0.05, which means ${\chi }_{\epsilon ,1}^2=3.84$.

The remainder of this section is organized as follows. First, we provide an evaluation of the accuracy achieved by the proposed approach (in comparison to other sub-pixel mapping approaches) using synthetic data. Then we provide an assessment and comparison using real hyperspectral scenes. The last section concludes with an evaluation of the impact of parameter settings on the newly developed JSSM approach.

4.1. Synthetic Experiments

4.1.1. Synthetic Image 1: HJ-1A

To simulate a real world situation under a fully controlled scenario, we used a real hyperspectral data set which was obtained by a Chinese environmental satellite, HJ-1A [55]. The HJ-1A satellite has a hyperspectral sensor with 115 spectral bands (in the spectral range from 0.45 to 0.95 μm) and spatial resolution of 100 m per pixel. The utilized HJ-1A image (150 × 150 pixels) was acquired on 19 August 2009 and was used as the original image in our experiments. The study site is located in Hanchuan City, Hubei Province, central China, and its surrounding area. Four land-cover classes, i.e., urban, agricultural land, water, and vegetation are used to characterize this image as Figure 3b shows. This image was degraded artificially to obtain a low-resolution synthetic image given a scale factor of three in this experiment.

As illustrated in Figure 3, the abundance map used as initial condition can severely impact the final sub-pixel mapping results. Because of this, it can be observed for the FCLS based sub-pixel mapping results that there are small pepper and salt patches existing for the whole image. For the HNN, MRF, and JSSM in which the abundance fractions can be altered, better results can be obtained. It still can be observed that many algorithms can generate similar results. Like results of HNN and MRF, our method can also eliminate some details due to model error and excessive smoothing. There is a trade-off between smoothing and detail preservation, in the sense that smoothing can restrain the model error while potentially eliminating small features. However, this kind of trade-off is inevitable owing to the nature of sub-pixel mapping, which aims to generate a finer classification map from the low-resolution remotely sensed image.

Table 1. Sub-pixel mapping results for the HJ-1A dataset.

Accuracy Indexes		Methods
		FCLS			SUnSAL-TV			MRF	JSSM
		AM	HC	HNN	AM	HC	HNN
IA (%)	Urban area	60.14	64.66	67.61	61.33	73.05	73.75	77.14	75.36
	Agricultural land	86.90	93.45	95.84	79.81	89.79	93.72	93.80	94.49
	Water	71.89	61.35	75.30	66.51	65.92	71.08	72.05	73.20
	Vegetation	51.49	56.67	50.26	55.86	61.26	57.54	56.43	60.03
AA (%)		67.61	69.03	72.25	65.88	72.51	74.02	74.86	75.77
OA (%)		64.23	65.82	69.42	63.60	71.30	72.55	74.22	74.33
Kappa		0.504	0.523	0.571	0.493	0.589	0.607	0.628	0.632
RMSE		0.192			0.180			Null
McNemar’s Test		901.3	1068.4	336.2	962.7	236.4	55.65	0.187	Null

shows a quantitative comparison of the AM, HC, HNN, MRF and the proposed JSSM methods. It is obvious that the sub-pixel mapping results of traditional methods are greatly affected by the spectral unmixing errors. It should be noted that the relationship between the accuracies of spectral unmixing and sub-pixel mapping is not monotonic as expected by many researches if the abundance fractions acts as a constraint. Although SUnSAL-TV can generate a better abundance map (RMSE value is 0.180) than that of FCLS, the sub-pixel mapping results exhibit worst performance when the AM method is utilized. Therefore, it is reasonable to develop sub-pixel mapping algorithms, which are conducted on the original remotely sensed image, such as MRF and JSSM. It can be seen that the proposed JSSM shows an improvement when compared with the traditional sub-pixel mapping methods, regardless of which spectral unmixing methods are utilized to generate the abundance map. When compared with the MRF, a slight improvement can also be observed in all quantitative indexes. For example, a gain of 0.91% in AA is obtained for the proposed method from the MRF method. Moreover, the McNemar’s test is a useful tool for determining if two classification methods have significantly different prediction rates. From Table 1, it can be seen that most of the values of the McNemar’s test, except MRF, are greater than the critical value (3.84). This implies that the proposed method has significantly different prediction rates compared with most traditional algorithms. Generally, the proposed JSSM can exhibit better performance by considering the endmember variability in sub-pixel mapping for a remotely sensed image.

4.1.2. Synthetic Image 2: Flightline C1 (FLC1)

In this section, we use an aerial data set with agricultural crop species and land use (the band number is 12) by an optical mechanical line scanner referred to as the University of Michigan M-7 system. The flightline used in this experiment (called FLC1) was collected on 28 June 1966 [56]. It was taken over the southern part of Tippecanoe County of Indiana. The size in pixels of the image is 80 × 160 pixels and the low-resolution image was generated using a resize factor of four as Figure 4a shows. Then a classification result was obtained for FLC1 by classifying the original high-resolution image with the commercial eCognition software. A total of eight land cover classes can be distinguished in Figure 4b–h give the sub-pixel mapping results.

As illustrated in Figure 4, it is obvious that the impact of the abundance map is quite important for the subsequent sub-pixel mapping procedure for traditional methods. As most land classes are kinds of agricultural crop species, it is easier to observe higher between-class similarity and the within-class dispersion. Therefore, the sub-pixel mapping results would be more susceptible to endmember variability. It can be seen that many pixels are misclassified by the MRF in which the endmember variability is not considered and the proposed JSSM method performs better. Moreover, owing to the special characteristic of the FLC1 dataset that contains few tiny features, the JSSM method can generate an acceptable visual result by imposing spatial regularization.

Table 2. Sub-pixel mapping results for the FLC1 dataset.

Accuracy Indexes		Methods
		FCLS			SUnSAL-TV			MRF	JSSM
		AM	HC	HNN	AM	HC	HNN
IA (%)	Red Cover	76.53	86.01	88.49	81.55	94.63	93.97	95.39	96.81
	Oats	58.95	74.70	80.73	55.07	94.23	79.58	85.83	94.70
	Wheat	93.80	94.94	96.16	94.20	99.35	98.20	98.20	99.67
	Soybeans	78.63	92.14	92.53	77.09	91.48	93.46	93.68	96.65
	Hay	41.74	55.51	56.89	43.89	86.91	73.32	65.24	88.84
	Pasture	62.42	75.88	78.79	62.52	77.43	76.95	80.58	85.71
	Alfalfa	65.54	78.24	75.93	65.54	79.66	74.78	83.30	82.42
	Corn	87.60	97.24	94.49	94.88	97.24	95.28	98.62	98.23
AA (%)		70.65	81.83	83.00	71.84	90.12	85.69	87.61	92.88
OA (%)		67.14	78.99	80.84	67.83	89.27	84.27	85.54	92.39
Kappa		0.623	0.757	0.779	0.631	0.876	0.818	0.833	0.912
RMSE		0.142			0.134			Null
McNemar’s Test		2893.4	1463.8	1193.1	2742.5	216.7	736.4	513.3	Null

shows a quantitative comparison of the AM, HC, HNN, MRF and the proposed JSSM methods. Different from the HJ-1A dataset mentioned above, it should be noted that SUnSAL-TV based sub-pixel mapping results show a better performance than that of FCLS here with a lower RMSE for the abundance map. From the quantitative results of HJ-1A and FLC1 datasets, the conclusion can therefore be drawn that the abundance map with higher accuracy cannot assure a better sub-pixel mapping result if the abundance fractions acts as the constraint. By incorporating this rich information in the original remotely sensed image, a great improvement in quantitative statistics can be observed for JSSM. The quantitative indexes lead to the same conclusion as the visual assessment for JSSM compared with MRF and the proposed JSSM exhibits better accuracies for AA, OA, and Kappa indexes. The OA increases from 85.54% to 92.39% for the JSSM compared with the MRF which results from the consideration of endmember variability. The McNemar’s test also demonstrates that JSSM gives significant different results compared to the critical value (3.84). Therefore, it can be concluded that the proposed JSSM performs better than other sub-pixel mapping methods from visual and quantitative assessments.

4.1.3. Synthetic Image 3: AVIRIS Indian Pines

The third image used in experiments is the commonly used AVIRIS Indian Pines dataset. Ten land cover classes were considered for classification. The utilized image is composed of 144 × 144 pixels and 220 bands, the ground truth data is also provided as the reference. The scale factor was set as three and Figure 5a shows the degraded AVIRIS hyperspectral image cube; Figure 5b shows the ground truth data in which 10 major land-cover classes can be distinguished and Figure 5c–h illustrate the sub-pixel mapping results with different algorithms.

The utilized Indian Pines image is well known for its land cover class spectral variability and severe noises in fields such as image classification and denoising. Therefore, the utilization of the endmember library can play an important role in the sub-pixel mapping procedure. For the Indian Pines image, the same conclusions can be drawn as for the previous two datasets. For the abundance map based sub-pixel mapping results, many isolated small false features can be found owing to the errors in the abundance map. While for the JSSM methods, these kinds of errors can be mitigated. Compared with the reference map, it can be seen that the proposed JSSM method can exhibit better performance by incorporating various endmember spectra of land cover classes than that of MRF. As shown in Figure 5, it can be distinguished that many patches of class Soybeans-clean are badly contaminated by other classes as the result of MRF while in the JSSM result, more homogeneous areas can be observed.

The sub-pixel mapping accuracies obtained in this experiment are listed in

Table 3. Sub-pixel mapping results for the Indian Pines dataset.

Accuracy Indexes		Methods
		FCLS			SUnSAL-TV			MRF	JSSM
		AM	HC	HNN	AM	HC	HNN
IA (%)	Corn-notill	44.00	48.73	49.92	65.46	79.98	83.70	53.72	89.86
	Corn-min	12.52	15.38	5.43	55.05	66.97	66.82	6.18	79.03
	Grass/Pasture	6.79	4.96	2.61	86.42	92.69	89.56	2.61	95.04
	Grass/Trees	52.87	53.22	73.91	82.09	92.00	98.96	85.91	96.00
	Hay-windrowed	92.70	99.16	100	94.38	100	100	100	100
	Soybeans-notill	61.86	66.27	64.80	56.55	83.36	84.68	66.57	80.71
	Soybeans-min	25.58	33.02	35.60	49.89	70.70	75.19	39.70	67.85
	Soybeans-clean	39.56	46.44	41.28	64.13	87.47	92.38	46.68	98.28
	woods	94.78	100	99.90	86.07	100	100	99.81	100
	Bldg	48.68	56.23	61.89	52.08	61.13	63.77	73.58	63.40
AA (%)		47.93	52.34	53.53	69.21	83.43	85.51	57.48	87.02
OA (%)		46.41	51.30	52.54	65.97	81.93	84.50	56.04	84.77
Kappa		0.391	0.443	0.457	0.611	0.792	0.821	0.496	0.825
RMSE		0.213			0.127			Null
McNemar’s Test		2486.6	2180.2	2026.1	965.9	77.9	0.52	1727.7	Null

to evaluate the effectiveness of the proposed method in quantitative fashion. It can be seen that the SUnSAL-TV based results show a great improvement compared with FCLS based results for both the abundance map and sub-pixel maps owing to the consideration of class spectral variability. The SUnSAL-TV based methods can even exhibit comparable results to the JSSM method. Moreover, it is easy to notice the big gaps between the JSSM and MRF methods for all quantitative indexes. For example, the JSSM improves the Kappa coefficient from 0.496 to 0.825 when compared with MRF. Generally, a better result can be obtained with the proposed JSSM quantitatively, especially for those images suffering from severe conditions of endmember variability.

4.2. Real Experiment-Nuance Dataset

To evaluate the practical application of the proposed method, a real experiment was implemented by acquiring a real hyperspectral image and a higher resolution color image for the same area, simultaneously. The original (low resolution) hyperspectral image (80 × 80 pixels) was collected using the Nuance NIR imaging spectrometer. The acquired hyperspectral image has 46 bands, and the spectral range is from 650–1100 nm, with 10 nm spectral sampling interval. The higher resolution color image (160 × 160 pixels) was obtained by a digital camera for the same scene, and the considered scale factor was two. The reference classification map was obtained by classifying the high resolution color image using the SVM. Four major land-cover classes can be distinguished in this experiment: soil, fresh vegetation, withered vegetation, and white paper. Figure 6a–c illustrates the original hyperspectral image, the high resolution color image, and the high resolution classification map used as reference, respectively. Figure 6d–i shows the sub-pixel mapping results obtained using AM, HC, HNN, MRF, and the proposed JSSM approach, respectively.

As shown in Figure 6, it is obvious that the proposed method can provide a better visual result compared with FCLS based sub-pixel mapping algorithms. Moreover, it can be seen that the difference between the JSSM and some other results is not so distinct. It is because the utilized Nuance image was acquired with a ground based imaging spectrometer and therefore, the impact of atmospheric and illumination conditions on the spectral variability can be suppressed. However, it still can be distinguished from the Figure 6 that the proposed JSSM shows a better result in some regions such as the shape of the paper at the bottom left of this image.

The accuracies and statistics of the different methods are listed in

Table 4. Sub-pixel mapping results for the Nuance dataset.

Accuracy Indexes		Methods
		FCLS			SUnSAL-TV			MRF	JSSM
		AM	HC	HNN	AM	HC	HNN
IA (%)	Soil	72.50	67.77	84.96	69.20	84.19	84.69	85.36	79.62
	Fresh vegetation	87.94	96.97	94.58	83.58	95.65	95.31	94.94	96.61
	Withered vegetation	70.03	78.52	77.73	72.84	80.28	79.33	80.66	86.01
	Paper	78.23	86.31	84.03	75.08	86.03	86.33	85.60	87.73
AA (%)		77.18	82.39	85.33	75.18	86.54	86.42	86.64	87.49
OA (%)		77.79	82.53	86.20	75.52	87.23	87.12	87.36	87.75
Kappa		0.699	0.764	0.812	0.669	0.827	0.825	0.828	0.834
RMSE		0.190			0.160			Null
McNemar’s Test		1508.9	808.7	84.4	1927.5	17.3	17.2	5.85	Null

. As illustrated in this table, the proposed JSSM also provides improvements in this experiment in terms of all indexes when compared with other sub-pixel mapping results. Although the MRF method can generate a comparable result with JSSM owing to the inconspicuous endmember variability for the Nuance dataset, the JSSM method still increases the AA from 86.64% to 87.49%, a gain of 0.85%, when compared with the MRF method. Moreover, the McNemar’s value indicates that JSSM gives a significantly different performance, compared to the other methods. Generally, both the visual assessment and quantitative accuracies indicate that JSSM performs better than the other sub-pixel mapping methods tested in this work.

4.3. Discussion

In the proposed JSSM method, three parameters were taken into consideration: the spatial constraint parameter ${\lambda }_{TV}$, the sparsity constraint parameter $\lambda$ and the penalty parameter $\mu$. These parameters have a great impact on the objective criteria’s optimization, or are sensitive to different data and significantly influence the sub-pixel mapping accuracy. To ensure our conclusions can be generalized to datasets with different gray level ranges, all images considered in the experiments were normalized to the range 0–1. We let all parameters vary from 10⁻⁶ to 10⁻¹ at a rate of 5 and Figure 7 and Figure 8 depict the impact of different values of parameters on the obtained accuracy for all considered experiments. Besides, the performance of the proposed JSSM method under different scale factors were also tested.

4.3.1. The Impact of the Penalty Parameter $\mu$

To evaluate the impact of the penalty parameter $\mu$, the highest accuracy of different $\mu$ values is obtained for all experimental datasets as Figure 7 shows. It can be seen that the trends for different datasets are not so consistent with each other. For the simulated FLC1 and Indian Pines images, the accuracy will decrease greatly since a certain threshold is satisfied. While for HJ-1A, the curve is more flat and a higher accuracy is further achieved with a greater value of $\mu$ for the Nuance dataset. However, acceptable results can be generated with a relative small value of $\mu$ for all datasets as shown in Figure 7. Therefore, the selection of the penalty parameter can be limited to around 10⁻³ empirically.

4.3.2. The Impact of the Regularization Parameter ${\lambda }_{TV}$ and $\lambda$

For the two regularization parameters, the optimal penalty parameter is determined firstly for different datasets. Figure 8 illustrates the accuracies of different parameters ${\lambda }_{TV}$ and $\lambda$ with all datasets. A similar pattern can be observed for simulated experiments and some curves are overlapped. It is easy to realize the difference of curves with simulated and real datasets; however, the same conclusions can be obtained that relative small values of ${\lambda }_{TV}$ and $\lambda$ can also give acceptable sub-pixel mapping results, while for the Indian Pines and Nuance datasets, a greater value of $\lambda$ can make the accuracy decrease dramatically. Generally, the values of ${\lambda }_{TV}$ and $\lambda$ can be set around 10⁻⁵ to ensure the generalization.

4.3.3. The Impact of Different Scale Factors

In the previous experiments, the scale factor was set as two, three, or four for different datasets. To illustrate the performance of the proposed JSSM under different scale factors, the simulated Indian Pines dataset was used to test different sub-pixel mapping algorithms. Figure 9 shows the accuracies of different methods when the scale factor is set as two, three, four, and six respectively. It is obvious that the proposed JSSM method generates the highest accuracies for different scale factors compared with other sub-pixel mapping methods. Moreover, it can also be found that the relationship between the accuracy of sub-pixel mapping and the scale factor is not monotonic. The sub-pixel mapping algorithms always achieve the highest accuracies when a moderate scale factor is given, such as three or four.

5. Conclusions and Future Lines

In this paper, we presented a new joint sparse sub-pixel mapping (JSSM) model for remotely sensed imagery. The proposed approach incorporates endmember variability into the sub-pixel mapping process by selecting variable representative spectra for each endmember in a library. Although the sub-pixel mapping is commonly regarded as the subsequent procedure of spectral unmixing, it is also shown in this paper that the use of a highly accurate abundance map does not always ensure a better sub-pixel mapping result. As a result, it is very important to address the problem of endmember variability for sub-pixel mapping on the original remotely sensed imagery directly, as opposed to traditional sub-pixel mapping algorithms which act on the abundance map. Compared with these traditional algorithms, which utilize a single spectrum for each endmember, the newly proposed JSSM can mitigate the impact of spectral variability and therefore provide a better sub-pixel mapping result. Our experimental results, conducted with a comprehensive database of synthetic and real images, indicate that JSSM is indeed an efficient sub-pixel mapping technique. Moreover, the proposed JSSM was also compared with its simplified version in which the endmember variability was not considered and the same conclusion can be drawn that the utilization of endmember variability improves the performance of sub-pixel mapping greatly. Although the proposed JSSM has been observed to outperform other traditional sub-pixel mapping methods, this model should be developed to incorporate additional datasets owing to the ill-posed nature of sub-pixel mapping. In this way, more robust and reliable sub-pixel mapping can be further generated. Additionally, it can also be seen in the experimental results that JSSM will give stairs (piecewise constant regions) owing to the utilization of the TV model. Our future work will focus on the extension of the proposed JSSM method for multi-temporal remotely sensed images in which the spectral variability is a severe problem owing to environmental, atmospheric, and temporal factors. Moreover, the automatic generation of the endmember library should be further taken into consideration and other spatial models, such as the non-local TV, will also be tested.