Nonlocal Total Variation Subpixel Mapping for Hyperspectral Remote Sensing Imagery

Ruyi Feng 1,2, Yanfei Zhong 1,2,*, Yunyun Wu 3, Da He 1,2, Xiong Xu 4 and Liangpei Zhang 1,2 1 State Key Laboratory of Information Engineering in Surveying, Mapping, and Remote Sensing, Wuhan University, Wuhan 430079, China; fry1988@whu.edu.cn (R.F.); whhd1992@163.com (D.H.); zlp62@whu.edu.cn (L.Z.) 2 Collaborative Innovation Center of Geospatial Technology, Wuhan University, Wuhan 430079, China 3 The Second Surveying and Mapping of Zhejiang Province, Hangzhou 310012, China; wenwubei@whu.edu.cn 4 College of Surveying and Geo-Informatics, Tongji University, Shanghai 200092, China; xvxiong@tongji.edu.cn * Correspondence: zhongyanfei@whu.edu.cn; Tel.: +86-027-6877-9969


Introduction
Due to the impact of the sensors' instantaneous field of view (IFOV) and the diversity of land-cover objects, mixed pixels are common in hyperspectral remote sensing images [1][2][3].Spectral unmixing can effectively solve the mixed pixel problem by acquiring the spectra of the endmembers and the abundance fraction images of these endmembers [4][5][6][7].However, the spatial distribution of each endmember or land-cover class in the pixel remains unknown [8,9], which is quite important in real applications such as subpixel classification and subpixel target detection.To solve this problem, the subpixel mapping technique was proposed to determine the subpixel location of each class or endmember based on the fractional abundances [10][11][12][13][14][15].In recent years, subpixel mapping has been applied in many practical applications, including land-cover change detection [16], waterline mapping [17], and landscape pattern indices calculation [18].Subpixel mapping methods have been proposed to handle the mapping problem inside pixels.They aim to generate a finer classification result from a low-resolution abundance map by maximizing the spatial dependence, which is the basic principle of many subpixel mapping methods [19].Subpixel mapping algorithms predict the spatial distribution of the different classes by dividing pixels into smaller subpixels based on the abundance values of each endmember, and then assign each subpixel to a single land-cover class, which can greatly enhance the spatial resolution of images.A number of subpixel mapping methods have been proposed based on spatial dependence, which refers to the fact that observations that are spatially close tend to be more similar than those that are further apart.Many subpixel mapping methods, both traditional and state-of-the-art, belong to this category, such as the linear optimization technique [20], the spatial attraction model (SASM) [21], the pixel-swapping subpixel mapping algorithm (PSSM) [22], neural networks [23,24], random fields [13,25], the geometric subpixel mapping algorithm (GSM) [26], evolutionary computing based subpixel mapping [27][28][29][30], and multi-agent systems [31].
Since the subpixel mapping issue is a typical ill-posed inverse problem whose aim is to predict unknown detail information at the subpixel scale from coarser pixels, subpixel mapping algorithms based on a spatial regularization prior have recently been studied by many researchers [12,14,32,33], as it is standard to use a regularization technique to make the original ill-posed problem well-posed [34].These algorithms transform the ill-posed subpixel mapping problem into a well-posed regularization problem by assuming spatial prior information about the unknown subpixel mapping result based on a maximum a posteriori (MAP) model.A finer-resolution map can then be reconstructed from a series of lower-resolution abundance fraction images.Different spatial prior models have been applied for subpixel mapping, such as the Laplacian prior, the total variation prior, and the bilateral total variation prior.These spatial prior models, such as the subpixel mapping based on total variation (TVSM) model proposed in [12], incorporate the local spatial information in the subpixel mapping model, and they have achieved better subpixel mapping performances.However, both the total variation regularization method and the Laplacian prior model just consider the local spatial homogeneity of the neighborhood system (horizontal and vertical) and reconstruct the main geometrical configuration, but they do not utilize all the potential spatial information [35], especially in the nonlocal spatial neighborhood.
In this paper, a new spatial regularization subpixel mapping algorithm based on a nonlocal prior model, namely nonlocal total variation subpixel mapping (NLTVSM), is proposed.In NLTVSM, a nonlocal total variation prior model [36] is designed as the spatial regularizer by incorporating the nonlocal spatial information in the subpixel mapping model.The nonlocal means-based regularization has been applied in many real applications, including hyperspectral image denoising [36] and noise estimation [37], super-resolution reconstruction [38], image segmentation [39], image classification [40], sparse unmixing [41], and so on.Differing from the previous TVSM method, which utilizes the spatial information in a limited local window, the proposed NLTVSM method combines the nonlocal spatial information of the whole image and the total variation framework.NLTVSM can use all possible self-predictions in the image to take advantage of the high degree of redundancy of the image, and a variational framework-based nonlocal operator is utilized as the spatial prior by averaging all the pixels in a sliding window and then computing the problem as a weighted graph gradient [42].Compared with the previous subpixel mapping algorithms, the experimental results using a simulated image, two synthetic hyperspectral remote sensing images, and a real remote sensing image demonstrate that the proposed NLTVSM algorithm can obtain an improved subpixel mapping accuracy, which is a little better than the latest TVSM algorithm.
The rest of this paper is organized as follows.Section 2 presents the basic idea of the proposed NLTVSM model.Section 3 provides the experimental results and analyses with simulated, synthetic, and real hyperspectral images.The conclusion is drawn in Section 4.

The Nonlocal Total Variation Subpixel Mapping Model for Hyperspectral Imagery
Differing from the previous TVSM [12], in this paper, the nonlocal total variation spatial operator is introduced into the subpixel mapping model to predict the fine structure, details, and texture, to enhance the subpixel mapping results.
In this part, the framework of the spatial regularization subpixel mapping model is first introduced, and the flowchart of the proposed NLTVSM is shown in Section 2.1.The spatial regularization method, the nonlocal total variation spatial operator, is presented in detail in Section 2.2.Finally, the proposed spatial regularization subpixel model, the NLTVSM algorithm, is described in Section 2.3.

Spatial Regularization Subpixel Mapping
Subpixel mapping techniques aim to obtain a higher-resolution classification image based on the low-resolution abundances, and divide a mixed pixel into several subpixels with a single land-cover class.The subpixel mapping problem involves estimating the possible distribution of every land-cover class inside a pixel at a subpixel scale, and the solution should not be unique, which is a typical ill-posed inverse problem.A simple example of subpixel mapping with three classes is displayed in Figure 1.In this example, the scale factor is 5, which means that one pixel in the original low-resolution image is divided into 5 ˆ5 subpixels.Figure 1a shows the original low-resolution fraction image, and the red, green, and blue pixels represent pure land-cover classes, while the white ones represent mixed pixels.The aim of subpixel mapping is to estimate the optimal distributions at the subpixel scale.Figure 1b,c show two possible distributions.Based on the spatial dependence rule, Figure 1b is more reasonable than the latter.
Remote Sens. 2016, 8, 250 3 of 20 regularization method, the nonlocal total variation spatial operator, is presented in detail in Section 2.2.Finally, the proposed spatial regularization subpixel model, the NLTVSM algorithm, is described in Section 2.3.

Spatial Regularization Subpixel Mapping
Subpixel mapping techniques aim to obtain a higher-resolution classification image based on the low-resolution abundances, and divide a mixed pixel into several subpixels with a single land-cover class.The subpixel mapping problem involves estimating the possible distribution of every land-cover class inside a pixel at a subpixel scale, and the solution should not be unique, which is a typical ill-posed inverse problem.A simple example of subpixel mapping with three classes is displayed in Figure 1.In this example, the scale factor is 5, which means that one pixel in the original low-resolution image is divided into 5 × 5 subpixels.Figure 1a shows the original low-resolution fraction image, and the red, green, and blue pixels represent pure land-cover classes, while the white ones represent mixed pixels.The aim of subpixel mapping is to estimate the optimal distributions at the subpixel scale.To better solve the ill-posed problem, a regularization technique can be considered to make it well-posed.In this way, the subpixel mapping problem can be transformed into a regularization problem, and a MAP model is utilized to regularize the subpixel mapping problem to be well-posed.The basic spatial regularization subpixel mapping model can be formulated as the following minimization problem: where ∈ × denotes the input low-resolution abundance images, with N pixels with n dimensions, where n is the number of land-cover classes.∈ ( × × ) × is the high-resolution subpixel mapping distribution map with s*s*N subpixels, and s is the scale factor.In the process of subpixel mapping, each pixel in the low-resolution abundance images corresponds to specific s*s subpixels.∈ ×( × × ) is the downsampling matrix, and the elements' values of matrix D are 0 and 1/(s 2 ).
In Equation ( 1), the first term describes the data fidelity term, and the second term is the spatial prior term.U(X) represents the spatial regularization term, and the parameter λ is the spatial regularization parameter which controls the balance of the data fidelity and regularization terms.Different spatial regularization models have previously been successfully applied for subpixel mapping, such as the Laplacian prior model, the total variation model, and the bilateral total variation (BTV) method.In this paper, the nonlocal total variation model is utilized to consider the To better solve the ill-posed problem, a regularization technique can be considered to make it well-posed.In this way, the subpixel mapping problem can be transformed into a regularization problem, and a MAP model is utilized to regularize the subpixel mapping problem to be well-posed.The basic spatial regularization subpixel mapping model can be formulated as the following minimization problem: where Y P R Nˆn denotes the input low-resolution abundance images, with N pixels with n dimensions, where n is the number of land-cover classes.X P R ps ˆs ˆNq ˆn is the high-resolution subpixel mapping distribution map with s*s*N subpixels, and s is the scale factor.In the process of subpixel mapping, each pixel in the low-resolution abundance images corresponds to specific s*s subpixels.D P R Nˆps ˆs ˆNq is the downsampling matrix, and the elements' values of matrix D are 0 and 1/(s 2 ).
In Equation (1), the first term describes the data fidelity term, and the second term is the spatial prior term.U(X) represents the spatial regularization term, and the parameter λ is the spatial regularization parameter which controls the balance of the data fidelity and regularization terms.Different spatial regularization models have previously been successfully applied for subpixel mapping, such as the Laplacian prior model, the total variation model, and the bilateral total variation (BTV) method.In this paper, the nonlocal total variation model is utilized to consider the spatial information in subpixel mapping, and the flowchart of the proposed NLTVSM method is shown in Figure 2. In the proposed method, the important improvement is the application of the nonlocal total variation method for the spatial consideration in the process of subpixel mapping.To better understand the proposed algorithm, the nonlocal total variation model should first be introduced.

The Nonlocal Total Variation Model
Unlike TVSM, in the proposed NLTVSM algorithm, the nonlocal total variation model is utilized as a spatial regularization term, which deals with the unknown pixels by searching for the similar pixels in a searching window to replace the current unknown pixels, with the aim being to exploit the similar patterns or patches in the image.The nonlocal-based method averages the other pixels with patches similar to that of the current one, as shown in Equation (2) [36]: where X describes the original image, x is the current pixel in the image, y is the central pixel in the similar window, ( ) * | ( • ) ( • )| ( ) ( ) is the normalizing factor, Gσ represents a Gaussian kernel with standard deviation σ, and h is a decay parameter., , and are the neighborhoods or the similar windows of pixels x, y, and z, respectively.
Figure 3 describes the basic principle of the nonlocal means model.In Figure 3, the whole image is assumed to be the searching window, which is utilized to restrict the search of similar windows in a real application.The small window centered at pixel P is the current window, and the windows centered at pixels q1, q2, and q3 are the similar windows to the current one.The nonlocal means operator is used to estimate the value of P using an average of all the pixels similar to the current one, such as the windows centered at pixels q1, q2, and q3, based on the weights of the measurement of the similarity w(P,q1), w(P,q2), and w(P, q3).In the proposed method, the important improvement is the application of the nonlocal total variation method for the spatial consideration in the process of subpixel mapping.To better understand the proposed algorithm, the nonlocal total variation model should first be introduced.

The Nonlocal Total Variation Model
Unlike TVSM, in the proposed NLTVSM algorithm, the nonlocal total variation model is utilized as a spatial regularization term, which deals with the unknown pixels by searching for the similar pixels in a searching window to replace the current unknown pixels, with the aim being to exploit the similar patterns or patches in the image.The nonlocal-based method averages the other pixels with patches similar to that of the current one, as shown in Equation (2) [36]: where X describes the original image, x is the current pixel in the image, y is the central pixel in the similar window, C pxq " e ´pGσ˚|Xpx `q´Xpz `q| 2 qpoq h 2 d pzq is the normalizing factor, G σ represents a Gaussian kernel with standard deviation σ, and h is a decay parameter.x `, y `, and z `are the neighborhoods or the similar windows of pixels x, y, and z, respectively.
Figure 3 describes the basic principle of the nonlocal means model.In Figure 3, the whole image is assumed to be the searching window, which is utilized to restrict the search of similar windows in a real application.The small window centered at pixel P is the current window, and the windows centered at pixels q1, q2, and q3 are the similar windows to the current one.The nonlocal means operator is used to estimate the value of P using an average of all the pixels similar to the current one, such as the windows centered at pixels q1, q2, and q3, based on the weights of the measurement of the similarity w(P,q1), w(P,q2), and w(P, q3).
is assumed to be the searching window, which is utilized to restrict the search of similar windows in a real application.The small window centered at pixel P is the current window, and the windows centered at pixels q1, q2, and q3 are the similar windows to the current one.The nonlocal means operator is used to estimate the value of P using an average of all the pixels similar to the current one, such as the windows centered at pixels q1, q2, and q3, based on the weights of the measurement of the similarity w(P,q1), w(P,q2), and w(P, q3).The nonlocal total variation (NLTV) functional is defined as a gradient-based function, which concerns both the nonlocal structures and textures, shown as: where X i (i P [1,n]) denotes one category of the distribution matrix X, and ∇ w X i is a nonlocal gradient which is defined as the vector of all the partial differences ∇ w X i (x,¨) at pixel x.X i (x) and X i (y) are the values located at pixels x and y in the ith distribution map X i .The weight w(x,y) is used to measure the similarity between pixels x and y.In Equation ( 3), there are two notations that should be given in detail.Firstly, the definition of the nonlocal gradient ∇ w X i (x,¨), which is calculated as shown in Equation ( 4) and, secondly, how to obtain the weight w(x,y), which is shown in Equation ( 5):  (5) where X i (N x ) and X i (N y ) respectively denote the square neighborhood of pixels x and y with a fixed window size on the category subpixel image X i .The Euclidean distance between these vectors, X i (y) ´Xi (x), expresses the similarity of the two pixels, x and y.In Equation (5), Z(x) is the normalizing constant, which is defined as Equation ( 6): 2,σ represents the Gaussian weighted Euclidean distance; σ is the standard deviation of the Gaussian kernel (σ > 0); and h acts as the degree of filtering, which controls the decay of the exponential function.

The NLTVSM Algorithm
In the previous section, we described the nonlocal total variation spatial operator in detail.Based on the basic spatial regularization subpixel mapping model in Equation ( 1), together with the nonlocal total variation model acting as the spatial regularization term, the NLTVSM algorithm can be written as Equation (7), where J NLTV (X) denotes the nonlocal total variation defined in Equation (3): As shown in Equation ( 7), the NLTVSM model has been built and the nonlocal total variation is designed as the spatial prior regularization term.Based on the equations introduced in Section 2.2, especially for the nonlocal total variation spatial operator, the objective function of the NLTVSM algorithm can be formulated.In this section, the Bregmanized operator splitting (BOS) algorithm [34] is used to optimize the NLTVSM model.The basic idea of the optimization algorithm can be stated as follows.
Firstly, the original NLTVSM optimization problem is enforced with the Bregman iteration [43] process as follows: The operator splitting technique is then used to solve the unconstrained subproblem in Equation ( 8) as follows: for j ě 0, X k+1,0 = X k , where j denotes the infinite inner iteration time: where δ is a positive value that satisfies 0 < δ < (2/||D T D||).The minimization problem in Equation ( 9) can be solved by the split Bregman method [44], with which the subproblem of X can be split into the following: where µ is the scale of the penalty term ||d-∇ w X i -b k || 2 , and its value is usually inversely proportional to the value of λ¨δ.The nonlocal total variation operator, denoted as J NLTV (X) in Equation ( 9), is rewritten as The solution of Equation ( 10), (X k+1 ,d k+1 ), can be acquired by alternating between the following two subproblems in Equation ( 11): The subproblem for X k+1 can be solved with the Euler-Lagrange equation, as follows: which provides: Therefore, X k+1 can be solved by the Gauss-Seidel algorithm [44].Meanwhile, the d k+1 in Equation ( 11) can be solved using a shrinkage operator as: The shrinkage operator is defined as follows: shrinkpp, tq " sign ppq max t|p| ´t, 0u To decrease the influence of the spectral unmixing errors, the winner-takes-all class determination strategy [13] is then utilized to obtain the final subpixel mapping results.

Experiments and Analysis
Experiments were conducted to test the performance of the proposed NLTVSM algorithm with one simulated image, two synthetic images, and one real image.The proposed algorithm was compared with five subpixel mapping algorithms: subpixel mapping based on a spatial attraction mode (SASM) [21], the pixel-swapping subpixel mapping algorithm (PSSM) [22], the genetic algorithm (GASM) [27], the geometric subpixel mapping algorithm (GSM) [26], the latest adaptive subpixel mapping method based on a MAP model, and a winner-takes-all class determination strategy (AMCDSM) with a total variation prior model [12] (denoted as TVSM in this paper).In addition, all the algorithms used the same fraction images for the subpixel mapping, which were obtained by the fully constrained least squares (FCLS) [45] method and probabilistic support vector machine (P-SVM) [46].
Since the subpixel mapping result is the same as the classification result, to evaluate the performance of the different subpixel mapping algorithms, the producer's accuracy (PA) (for single classes), the overall accuracy (OA), and the Kappa coefficient were adopted, and were used to verify the performance of the different classification algorithms with the help of ENVI software [47].For the simulated hyperspectral images and the synthetic images, the low-resolution hyperspectral remote sensing images were first obtained by downsampling the original high-resolution images with a fixed scale.A spectral unmixing technique was then utilized to obtain the unmixing results-the fractional abundance images-as the input of the subpixel mapping.With the same fractional abundance images as the input, the different subpixel mapping algorithms obtained different performances and obtained different subpixel mapping results.The reference classification image used to verify the performance of the different subpixel mapping results was obtained by classification of the original high-resolution images by a hard classification method, such as a support vector machine (SVM) [48] or the minimum distance hard classification algorithm [49], which were implemented in ENVI software.For the real hyperspectral image, the low-resolution (LR) hyperspectral image was collected with a Nuance NIR imaging spectrometer, and the high-resolution (HR) image was taken by a digital camera for the same area at the same time.The spectral unmixing algorithm was then used with the LR images to obtain the fraction images, and the hard classification algorithm was used to obtain the ground truth from the HR data.More information about the experimental datasets is provided in the following section.

Experimental Design and Datasets
The first experimental dataset (simulated image) with 400 ˆ400 pixels and 50 bands, as shown in Figure 4a, was generated using the USGS spectral library following the methodology in [13].This dataset consists of four typical land-cover classes: water, tree, agricultural field, and impervious layer.Figure 4b displays the true subpixel mapping image, acting as the reference image.In addition, a spectral unmixing algorithm was used with the low-resolution hyperspectral images (LR), which were acquired via a 4 ˆ4 mean filter from the original high-resolution simulated hyperspectral image, to obtain the fractional abundance images, as shown in Figure 4c.The subpixel scale factor (s) was set to 4 in this dataset.
Two synthetic hyperspectral remote sensing images were also used to test the performance of the different subpixel mapping algorithms.The first image was the Washington DC Mall Hyperspectral Digital Imagery Collection Experiment (HYDICE) image, with 200 ˆ300 pixels and 191 bands (as shown in Figure 5a), which has four classes of land cover, i.e., water, grass, tree, and road.The other image with more classes was the HYDICE Urban image with 300 ˆ300 pixels and 187 bands, as shown in Figure 6a, which is displayed in false color.The area used in the experiment consisted of six major classes of roof1, tree, concrete road, roof/shadow, grass, and asphalt road.
in Figure 4a, was generated using the USGS spectral library following the methodology in [13].This dataset consists of four typical land-cover classes: water, tree, agricultural field, and impervious layer.Figure 4b displays the true subpixel mapping image, acting as the reference image.In addition, a spectral unmixing algorithm was used with the low-resolution hyperspectral images (LR), which were acquired via a 4 × 4 mean filter from the original high-resolution simulated hyperspectral image, to obtain the fractional abundance images, as shown in Figure 4c.The subpixel scale factor (s) was set to 4 in this dataset.The classification maps for the evaluation of these two synthetic datasets were obtained from the original hyperspectral images with the SVM hard classification algorithm, and the reference classification images for the two datasets are shown in Figures 5c and 6c.In addition, the results of the classification with the SVM algorithm were also tested with the ROIs selected manually, as shown in Figures 5b and 6b, and the accuracies were 97.12% [32] and 97.03% for the OA values, corresponding to the Washington DC Mall HYDICE image (shown in Table 1) and the HYDICE Urban image (shown in Table 2), respectively.These accuracies were also considered to be sufficient for the subpixel mapping experiments.To simulate the coarse-resolution hyperspectral remote sensing images, a 4 × 4 mean filter was utilized for the downsampling, and the fractional abundances were obtained with the spectral unmixing algorithm, as shown in Figures 5d and 6d, as inputs.In these two experiments, the subpixel mapping scale factor (s) was set to 4. The classification maps for the evaluation of these two synthetic datasets were obtained from the original hyperspectral images with the SVM hard classification algorithm, and the reference classification images for the two datasets are shown in Figures 5c and 6c.In addition, the results of the classification with the SVM algorithm were also tested with the ROIs selected manually, as shown in Figures 5b and 6b, and the accuracies were 97.12% [32] and 97.03% for the OA values, corresponding to the Washington DC Mall HYDICE image (shown in Table 1) and the HYDICE Urban image (shown in Table 2), respectively.These accuracies were also considered to be sufficient for the subpixel mapping experiments.To simulate the coarse-resolution hyperspectral remote sensing images, a 4 × 4 mean filter was utilized for the downsampling, and the fractional abundances were obtained with the spectral unmixing algorithm, as shown in Figures 5d and 6d, as inputs.In these two experiments, the subpixel mapping scale factor (s) was set to 4. The classification maps for the evaluation of these two synthetic datasets were obtained from the original hyperspectral images with the SVM hard classification algorithm, and the reference classification images for the two datasets are shown in Figures 5c and 6c.In addition, the results of the classification with the SVM algorithm were also tested with the ROIs selected manually, as shown in Figures 5b and 6b, and the accuracies were 97.12% [32] and 97.03% for the OA values, corresponding to the Washington DC Mall HYDICE image (shown in Table 1) and the HYDICE Urban image (shown in Table 2), respectively.These accuracies were also considered to be sufficient for the subpixel mapping experiments.To simulate the coarse-resolution hyperspectral remote sensing images, a 4 ˆ4 mean filter was utilized for the downsampling, and the fractional abundances were obtained with the spectral unmixing algorithm, as shown in Figures 5d and 6d, as inputs.In these two experiments, the subpixel mapping scale factor (s) was set to 4. The real experimental dataset was acquired by a Nuance NIR imaging spectrometer, as shown in Figure 7a, and consisted of 46 bands, with 50 × 50 pixels.The spectral unmixing algorithm was applied to undertake the unmixing and obtain the fractional abundances, as shown in Figure 7b.To assess the performance of the different subpixel mapping algorithms, the HR image (Figure 7c)-150 × 150 pixels and three bands-was also collected for the same area using a digital camera at the same time.To verify the subpixel mapping results, the reference classification map was then obtained by the SVM algorithm, with the ROIs selected manually, as shown in Figure 7d,e, respectively.Furthermore, the accuracy of the SVM classification method for the Nuance image was also determined, as shown in Table 3.The OA was 99.27%, which confirmed that the result was suitable for testing the other subpixel mapping methods.In addition, the subpixel mapping scale factor (s) was set to three for the Nuance dataset.withered vegetable fresh vegetable background The real experimental dataset was acquired by a Nuance NIR imaging spectrometer, as shown in Figure 7a, and consisted of 46 bands, with 50 ˆ50 pixels.The spectral unmixing algorithm was applied to undertake the unmixing and obtain the fractional abundances, as shown in Figure 7b.To assess the performance of the different subpixel mapping algorithms, the HR image (Figure 7c)-150 ˆ150 pixels and three bands-was also collected for the same area using a digital camera at the same time.To verify the subpixel mapping results, the reference classification map was then obtained by the SVM algorithm, with the ROIs selected manually, as shown in Figure 7d,e, respectively.Furthermore, the accuracy of the SVM classification method for the Nuance image was also determined, as shown in Table 3.The OA was 99.27%, which confirmed that the result was suitable for testing the other subpixel mapping methods.In addition, the subpixel mapping scale factor (s) was set to three for the Nuance dataset.(e) (f) (g) (e) (f) (g) The subpixel mapping results for the simulated image are shown in Figure 8b-g.For convenience, the true subpixel mapping image is also shown in Figure 8a.It can be observed that the results obtained by SASM, PSSM, GASM, and TVSM display similar distributions to the true subpixel mapping image, but the result of GSM exhibits a great difference, especially for the simulated impervious layer of the linear feature in red.However, the NLTVSM algorithm performs the best in this simulated dataset, especially in the transverse line in the top of the image, wiping off the wrong mapping at the subpixel scale.Since the nonlocal total variation spatial operator can obtain more spatial prior information according to the larger nonlocal searching window, the proposed NLTVSM achieves an improved performance.

Experimental Results and Analysis
To evaluate the different performances quantitatively, Table 4 provides the quantitative assessment results of the above subpixel mapping methods for the simulated image.It can be seen that the OA and Kappa are the same for SASM, PSSM, and GASM.The performance of GSM is not as good, at only 89.04% for the OA and 0.84 for Kappa.For the spatial regularization subpixel mapping algorithms-TVSM and NLTVSM-they both obtain competitive accuracies.Owing to the nonlocal searching strategy considering more spatial information in the nonlocal neighborhood, which can find more similar spatial distribution patterns in a larger neighborhood, NLTVSM obtains the best quantitative performance.For the two synthetic hyperspectral images, the clear advantage of the proposed NLTVSM algorithm can be seen.For the results of the Washington DC Mall HYDICE image shown in Figure 9, it can be observed that the result of NLTVSM exhibits less salt-and-pepper noise, compared with the other methods, and the edge of the meadow and the lake line in the top-left corner show The subpixel mapping results for the simulated image are shown in Figure 8b-g.For convenience, the true subpixel mapping image is also shown in Figure 8a.It can be observed that the results obtained by SASM, PSSM, GASM, and TVSM display similar distributions to the true subpixel mapping image, but the result of GSM exhibits a great difference, especially for the simulated impervious layer of the linear feature in red.However, the NLTVSM algorithm performs the best in this simulated dataset, especially in the transverse line in the top of the image, wiping off the wrong mapping at the subpixel scale.Since the nonlocal total variation spatial operator can obtain more spatial prior information according to the larger nonlocal searching window, the proposed NLTVSM achieves an improved performance.
To evaluate the different performances quantitatively, Table 4 provides the quantitative assessment results of the above subpixel mapping methods for the simulated image.It can be seen that the OA and Kappa are the same for SASM, PSSM, and GASM.The performance of GSM is not as good, at only 89.04% for the OA and 0.84 for Kappa.For the spatial regularization subpixel mapping algorithms-TVSM and NLTVSM-they both obtain competitive accuracies.Owing to the nonlocal searching strategy considering more spatial information in the nonlocal neighborhood, which can find more similar spatial distribution patterns in a larger neighborhood, NLTVSM obtains the best quantitative performance.For the two synthetic hyperspectral images, the clear advantage of the proposed NLTVSM algorithm can be seen.For the results of the Washington DC Mall HYDICE image shown in Figure 9, it can be observed that the result of NLTVSM exhibits less salt-and-pepper noise, compared with the other methods, and the edge of the meadow and the lake line in the top-left corner show clear boundaries.
The quantitative results in Table 5 allow the same conclusion to be made.The proposed NLTVSM algorithm obtains a higher accuracy than the traditional methods and the state-of-the-art ones.For instance, the OA value of NLTVSM is 78.11% and the Kappa is 0.69 for the Washington DC Mall dataset.This is an improvement of 1.15% and 0.02 in the subpixel mapping accuracy (OA and Kappa coefficient, respectively) compared with the TVSM approach.The results for the HYDICE Urban image are presented in Figure 10b-g, and a detailed verification of the results is provided in Table 6.Here, it can be observed that all the algorithms display different performances.It can be seen that the flat zones, such as grass or asphalt road, consist of lots of scatter, perhaps due to the fraction errors derived from the spectral unmixing process.Compared with the SASM, PSSM, and GASM methods, the GSM, TVSM, and NLTVSM approaches exhibit smoother visual results, benefiting from the consideration of the spatial information.Table 6 shows consistent results.Compared with the SASM, PSSM, and GASM methods, GSM, TVSM, and NLTVSM possess obvious advantages, and show an improvement of at least 4.63% in the OA and 0.05 in Kappa.The NLTVSM algorithm obtains the best accuracy among these algorithms by a narrow margin: 70.12% for the OA and 0.62 for Kappa, an OA improvement of 2.09% compared with TVSM.The experimental results show that NLTVSM is better than TVSM for more complex images with more classes.In the real hyperspectral data experiment, the eight algorithms all obtain similar results for the withered vegetable class, as shown in Figure 11.The subpixel mapping results of the SASM, PSSM, and GASM algorithms with fraction constraints are poor, with significant salt-and-pepper noise.The other three subpixel mapping algorithms-GSM, TVSM, and the proposed NLTVSM-obtain better results with smoother and more continuous land-cover distributions.These algorithms relax the constraint of the fraction values, and can decrease the errors in the spectral unmixing.For the fresh vegetable and the background classes in particular, the subpixel mapping result of the NLTVSM method is better than the traditional algorithms.Comparing the proposed NLTVSM with the GSM and TVSM algorithms, although they all consider the spatial information as prior knowledge in the process of subpixel mapping, NLTVSM obtains a better result by utilizing the nonlocal total variation regularization operator, which can take advantage of the high degree of redundancy in the nonlocal spatial information of the image to predict more fine structure and details, and can suppress most of the salt-and-pepper noise, e.g., for the withered vegetable class.
The quantitative comparison using the OA and Kappa provided in Table 7 also allows the same observation to be made.As shown in Table 7, GSM, TVSM, and the proposed NLTVSM show a great improvement in subpixel mapping accuracy over the traditional subpixel mapping approaches.Compared with GSM and TVSM, NLTVSM obtains the highest subpixel mapping accuracies in Table 7.The OA and Kappa of NLTVSM are equal to 76.60% and 0.64%, respectively, for the real Nuance image.The main reason for the higher accuracy of the NLTVSM algorithm is as follows.TVSM only accounts for the spatial homogeneity of the first-order pixel neighborhood system.However, NLTVSM tries to make full use of the potential nonlocal spatial information by utilizing all possible self-similarities existing in the image.Based on the above analysis, the NLTVSM algorithm provides an effective option to perform the task of subpixel mapping for hyperspectral remote sensing imagery.

Sensitivity Analysis
(1) Discussion on sensitivity analysis for the spatial regularization parameter In the optimization problem of the NLTVSM algorithm, there is an important parameter, λ-the spatial regularization parameter-which controls the balance of the data fidelity and regularization terms, and significantly influences the final subpixel mapping results.Figure 12 analyzes the effect of this parameter for the Washington DC Mall and Nuance images.NLTVSM method is better than the traditional algorithms.Comparing the proposed NLTVSM with the GSM and TVSM algorithms, although they all consider the spatial information as prior knowledge in the process of subpixel mapping, NLTVSM obtains a better result by utilizing the nonlocal total variation regularization operator, which can take advantage of the high degree of redundancy in the nonlocal spatial information of the image to predict more fine structure and details, and can suppress most of the salt-and-pepper noise, e.g., for the withered vegetable class.The quantitative comparison using the OA and Kappa provided in Table 7 also allows the same observation to be made.As shown in Table 7, GSM, TVSM, and the proposed NLTVSM show a great improvement in subpixel mapping accuracy over the traditional subpixel mapping approaches.Compared with GSM and TVSM, NLTVSM obtains the highest subpixel mapping accuracies in Table 7.The OA and Kappa of NLTVSM are equal to 76.60% and 0.64%, respectively, for the real Nuance image.The main reason for the higher accuracy of the NLTVSM algorithm is as follows.TVSM only accounts for the spatial homogeneity of the first-order pixel neighborhood system.However, NLTVSM tries to make full use of the potential nonlocal spatial information by utilizing all possible self-similarities existing in the image.Based on the above analysis, the NLTVSM algorithm provides an effective option to perform the task of subpixel mapping for hyperspectral remote sensing imagery.

Sensitivity Analysis
(1) Discussion on sensitivity analysis for the spatial regularization parameter In the optimization problem of the NLTVSM algorithm, there is an important parameter, λ-the spatial regularization parameter-which controls the balance of the data fidelity and regularization terms, and significantly influences the final subpixel mapping results.Figure 12 analyzes the effect of this parameter for the Washington DC Mall and Nuance images.As shown in Figure 12, it can be seen that the best subpixel mapping result is acquired when the value of λ is equal to 0.05 and 0.5 for the two images, respectively.To achieve a good result and effectively exploit the useful nonlocal spatial information, appropriate regularization parameters should be found according to the different images or land-cover classes.
(2) The effect of different spectral unmixing approaches on NLTVSM In this paper, the fractional abundances, derived from the spectral unmixing, act as the input of the subpixel mapping.Since uncertainty exists for the different spectral unmixing methods, and different unmixing results lead to different subpixel mapping outputs, the sensitivity of the different spectral unmixing approaches is discussed in this part.As an example, the Washington DC Mall image is used for analyzing the sensitivity of the unmixing results for the subpixel mapping, and the spectral unmixing algorithms adopted for comparison are FCLS, FCLS, and P-SVM, and a sparse regression method, sparse unmixing via variable splitting, and augmented Lagrangian (SUnSAL) [50].Figure 13 shows the final subpixel mapping results.As shown in Figure 12, it can be seen that the best subpixel mapping result is acquired when the value of λ is equal to 0.05 and 0.5 for the two images, respectively.To achieve a good result and effectively exploit the useful nonlocal spatial information, appropriate regularization parameters should be found according to the different images or land-cover classes.
(2) The effect of different spectral unmixing approaches on NLTVSM In this paper, the fractional abundances, derived from the spectral unmixing, act as the input of the subpixel mapping.Since uncertainty exists for the different spectral unmixing methods, and different unmixing results lead to different subpixel mapping outputs, the sensitivity of the different spectral unmixing approaches is discussed in this part.As an example, the Washington DC Mall image is used for analyzing the sensitivity of the unmixing results for the subpixel mapping, and the spectral unmixing algorithms adopted for comparison are FCLS, FCLS, and P-SVM, and a sparse regression method, sparse unmixing via variable splitting, and augmented Lagrangian (SUnSAL) [50].Figure 13 shows the final subpixel mapping results.It can be seen that different spectral unmixing results lead to different subpixel mapping outputs.FCLS, which is one of the classical linear spectral unmixing approaches, and is widely used for its robustness and tractability, leads to a satisfactory subpixel mapping result.For the sparse unmixing method, SUnSAL obtains a similar subpixel mapping result, although the spectral unmixing result evaluated by the signal-to-reconstruction error (SRE) value [50] of SUnSAL with 5.15 dB is a little higher than FCLS with 4.96 dB.The possible reason for this is the spatial consideration in the nonlocal total variation subpixel mapping method, which results in the unmixing difference being reduced.The spectral unmixing methods adopted in this paper, FCLS and P-SVM, obtain the best subpixel mapping results.
(3) Discussion on the effect of adding a shade endmember In this paper, most of the final subpixel mapping results consist of a set of typical materials.For example, water, grass, tree, road, roof, etc. Shade is usually not considered as a typical endmember or land-cover component [51,52], but is treated as a variant of the brightness of other endmembers.However, it is interesting to discuss the effect of adding a shade endmember in the experiments, especially as water and shade have similar spectra.In this part, we discuss the effect Water Grass Tree Road It can be seen that different spectral unmixing results lead to different subpixel mapping outputs.FCLS, which is one of the classical linear spectral unmixing approaches, and is widely used for its robustness and tractability, leads to a satisfactory subpixel mapping result.For the sparse unmixing method, SUnSAL obtains a similar subpixel mapping result, although the spectral unmixing result evaluated by the signal-to-reconstruction error (SRE) value [50] of SUnSAL with 5.15 dB is a little higher than FCLS with 4.96 dB.The possible reason for this is the spatial consideration in the nonlocal total variation subpixel mapping method, which results in the unmixing difference being reduced.The spectral unmixing methods adopted in this paper, FCLS and P-SVM, obtain the best subpixel mapping results.
(3) Discussion on the effect of adding a shade endmember In this paper, most of the final subpixel mapping results consist of a set of typical materials.For example, water, grass, tree, road, roof, etc. Shade is usually not considered as a typical endmember or land-cover component [51,52], but is treated as a variant of the brightness of other endmembers.However, it is interesting to discuss the effect of adding a shade endmember in the experiments, especially as water and shade have similar spectra.In this part, we discuss the effect of adding a shade endmember to the Washington DC Mall image.The comparison results are shown in Figure 14.

Running Times of the Different Models for Different Images
The running times of the different models for the two synthetic images were tested on a PC with an Intel Core i3-2100 CPU @ 3.1 GHz and 8 GB RAM, using MATLAB R2011b.Since the different strategies used in the algorithms may lead to different running times, especially on the MATLAB platform, we just compared the algorithms proposed by our research group, i.e., TVSM and the proposed NLTVSM.All of the parameters and the iteration times were set to obtain the optimal results.The running times of the different models are listed in Table 9.It can be observed from Table 9 that the running time of the proposed method, NLTVSM, is a little less than the other algorithm, TVSM, owing to the fact that the core process of NLTVSM is coded in Visual C++ 6.0.Furthermore, the use of fast Fourier transform (FFT) enhances the computational efficiency of TVSM, and it also greatly improves the efficiency of TVSM.

Conclusions
In this paper, the nonlocal total variation subpixel mapping algorithm, namely NLTVSM, has been proposed for hyperspectral remote sensing imagery.In NLTVSM, a nonlocal total variation model is designed as the spatial regularization term to utilize the nonlocal spatial information.Differing from the previous spatial regularization subpixel mapping algorithm based on total variation (TVSM), NLTVSM extracts the nonlocal spatial information by the use of a variational framework based nonlocal operator averaging the set of subpixels in a certain size of sliding window in the image.This combines the nonlocal spatial information of the whole image and the total variation framework, and utilizes all possible self-predictions in the subpixel mapping images to take advantage of the high degree of redundancy of the image.Experimental results using a simulated image, two synthetic hyperspectral remote sensing images, and a real hyperspectral image demonstrate that the proposed NLTVSM performs better than the traditional subpixel mapping algorithms and the latest TVSM algorithm.NLTVSM can obtain a higher subpixel mapping precision, a more accurate distribution map with less salt-and-pepper effect at the subpixel scale, and is also efficient.Hence, it provides a new option to perform the task of subpixel mapping for hyperspectral remote sensing imagery.In our future work, an adaptive parameter selection method will be studied.

Figure 3 .
Figure 3.The basic principle of the nonlocal means model.

Figure 3 .Figure 3 .
Figure 3.The basic principle of the nonlocal means model.

Figure 4 .Figure 4 .Figure 5 .Figure 6 .
Figure 4. Simulated image.(a) Original simulated hyperspectral image (400 × 400); (b) reference classification image; and (c) fractional abundances (100 × 100).Two synthetic hyperspectral remote sensing images were also used to test the performance of the different subpixel mapping algorithms.The first image was the Washington DC Mall Hyperspectral Digital Imagery Collection Experiment (HYDICE) image, with 200 × 300 pixels and 191 bands (as shown in Figure5a), which has four classes of land cover, i.e., water, grass, tree, and road.The other image with more classes was the HYDICE Urban image with 300 × 300 pixels and

Figures 8 -
Figures 8-11 show the subpixel mapping results using SASM, PSSM, GASM, GSM, TVSM, and the proposed NLTVSM for the simulated image, the Washington DC Mall HYDICE image, the HYDICE Urban image, and the Nuance image, respectively.The visual comparisons in Figures 8-11 show the varying performances of the different subpixel mapping methods.

Figures 8 -Figure 8 .
Figures 8-11 show the subpixel mapping results using SASM, PSSM, GASM, GSM, TVSM, and the proposed NLTVSM for the simulated image, the Washington DC Mall HYDICE image, the HYDICE Urban image, and the Nuance image, respectively.The visual comparisons in Figures 8-11 show the varying performances of the different subpixel mapping methods.

Figure 12 .
Figure 12.Sensitivity analysis for the spatial regularization parameter λ.(a) Washington DC Mall image; and (b) real Nuance image.

Figure 12 .
Figure 12.Sensitivity analysis for the spatial regularization parameter λ.(a) Washington DC Mall image; and (b) real Nuance image.

Table 1 .
Accuracy of the classification method for the Washington DC Mall HYDICE image.

Table 2 .
Accuracy of the classification method for the HYDICE Urban image.

Table 2 .
Accuracy of the classification method for the HYDICE Urban image.

Table 3 .
Accuracy of the classification method for the Nuance dataset.

Table 3 .
Accuracy of the classification method for the Nuance dataset.

Table 4 .
Accuracy of the different algorithms with the simulated image.
Note: The highest accuracies in the table are marked with bold and underlined.

Table 4 .
Accuracy of the different algorithms with the simulated image.
Note: The highest accuracies in the table are marked with bold and underlined.

Table 5 .
Accuracy of the different algorithms with the Washington DC Mall image.
Note: The highest accuracies in the table are marked with bold and underlined.

Table 6 .
Accuracy of the different algorithms with the HYDICE Urban image.Note: The highest accuracies in the table are marked with bold and underlined.

Table 7 .
Accuracy of the different algorithms with the Nuance dataset.
Note: The highest accuracies in the table are marked with bold and underlined.

Table 7 .
Accuracy of the different algorithms with the Nuance dataset.
Note: The highest accuracies in the table are marked with bold and underlined.

Table 9 .
Running times of the different models for the different images.