Adaptive non-local Euclidean medians sparse unmixing for hyperspectral imagery

ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

Contents lists available at ScienceDirect

ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs

Adaptive non-local Euclidean medians sparse unmixing for hyperspectral imagery Ruyi Feng, Yanfei Zhong ⇑, Liangpei Zhang State Key Laboratory of Information Engineering in Surveying, Mapping, and Remote Sensing, Wuhan University, PR China

a r t i c l e

i n f o

Article history: Received 13 April 2014 Received in revised form 20 June 2014 Accepted 17 July 2014

Keywords: Hyperspectral imagery Non-local Euclidean medians Non-local means Adaptive Sparse unmixing

a b s t r a c t Sparse unmixing models based on sparse representation theory and a sparse regression model have been successfully applied to hyperspectral remote sensing image unmixing. To better utilize the abundant spatial information and improve the unmixing accuracy, spatial sparse unmixing methods such as the non-local sparse unmixing (NLSU) approach have been proposed. Although the NLSU method utilizes non-local spatial information as the spatial regularization term and obtains a satisfactory unmixing accuracy, the final abundances are affected by the non-local neighborhoods and drift away from the true abundance values when the observed hyperspectral images have high noise levels. Furthermore, NLSU contains two regularization parameters which need to be appropriately set in real applications, which is a difficult task and often has a high computational cost. To solve these problems, an adaptive non-local Euclidean medians sparse unmixing (ANLEMSU) method is proposed to improve NLSU by replacing the non-local means total variation spatial consideration with the non-local Euclidean medians filtering approach. In addition, ANLEMSU utilizes a joint maximum a posteriori (JMAP) strategy to acquire the relationships between the regularization parameters and the estimated abundances, and achieves the fractional abundances adaptively, without the need to set the two regularization parameters manually. The experimental results using both simulated data and real hyperspectral images indicate that ANLEMSU outperforms the previous sparse unmixing algorithms and, hence, provides an effective option for the unmixing of hyperspectral remote sensing imagery. Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.

1. Introduction Spectral unmixing was first proposed to cope with the problem of mixed pixels encountered in remote sensing imagery (Demarchi et al., 2012; Tits et al., 2012) and to identify the components of the mixed spectra (also called endmembers) in each pixel, together with their proportions (known as abundances) (Bioucas-Dias et al., 2012, 2013; Keshava and Mustard, 2002; Ma et al., 2014; Tong et al., 2014; Li et al., 2014). Due to the computational tractability and flexibility of the linear mixture model (LMM) (Pu et al., 2014), the LMM has been widely studied. All the unmixing models discussed in this paper are based on the LMM. The conventional spectral unmixing techniques usually consist of an endmember identification step (Zhong et al., 2014b) and an abundance inversion process (Qu et al., 2014). However, the identification of the endmembers is a challenging task because of the

⇑ Corresponding author. Tel.: +86 27 68779969. E-mail address: [email protected] (Y. Zhong).

lack of prior knowledge about the endmember signatures (Ma et al., 2014). Sparse unmixing is a novel spectral unmixing technique that is also known as the dictionary-based semi-blind hyperspectral unmixing approach, and it has been widely studied (Iordache, 2011; Iordache et al., 2011, 2012, 2014a,b; Wu et al., 2011; Zhao et al., 2013; Zhong et al., 2014a; Liu and Zhang, 2014; Zhu et al., 2014; Tang et al., 2014). Sparse unmixing assumes that the observed hyperspectral remote sensing imagery can be expressed in a linear sparse regression (Plaza et al., 2011). Unmixing via sparse representation can be reformulated as finding the best combination of endmembers in a prior large standard spectral library, and the sparsity constraint is imposed on the corresponding endmembers’ fractional abundances, since the number of endmembers in each mixed pixel is sparse in reality (Iordache et al., 2011). To solve this linear sparse spectral unmixing problem, the sparse unmixing via variable splitting and augmented Lagrangian (SUnSAL) algorithm (Iordache et al., 2011) was proposed. The application of a spectral library in the SUnSAL algorithm successfully circumvents the determination of the endmembers, and better

http://dx.doi.org/10.1016/j.isprsjprs.2014.07.009 0924-2716/Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.

10

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

results are achieved than with the traditional linear spectral unmixing methods such as non-negative constrained least squares (Iordache, 2011; Iordache et al., 2011). As the research into sparse unmixing has progressed, the spatial information existing in the remote sensing imagery has also been utilized, and several spatial sparse unmixing algorithms have been developed, such as sparse unmixing via variable splitting augmented Lagrangian and total variation (SUnSAL-TV) (Iordache et al., 2012) and non-local sparse unmixing (NLSU) (Zhong et al., 2014a). Compared with SUnSAL-TV, which utilizes the spatial information in a first-order pixel neighborhood system, NLSU has been proved to be an efficient and powerful spatial sparse unmixing algorithm that considers the non-local spatial information of the whole abundance image. In the NLSU algorithm, each mixed pixel’s abundance is reformulated with the non-local means (NLM) (Buades et al., 2005a,b) total variational method as the spatial consideration, using a weighted averaging of a group of nonlocal pixels’ abundances, and the Euclidean distance is utilized to measure the similarity of each neighborhood. However, due to the existence of the inaccurate estimated unmixing abundances, the weights measured by the Euclidean distance will not be exact. Unfortunately, the performance of NLSU largely depends on the reliability of the weights, and it is inevitable that weights computed from estimated fractional abundances will lead to an inaccurate spatial relationship. That is to say, the outliers in the weights, together with the corresponding unreliable similar windows, should be excluded in the process of computing the fractional abundances for the current pixel. Moreover, the two regularization parameters in the NLSU model, one before the sparsity regularization term, and the other connected with the spatial constraint for the contribution of the spatial regularization term, need to be determined in real applications, which is a difficult task and is quite time-consuming. To obtain a more reliable spatial correlation and to obtain the regularization parameters of the model adaptively, an adaptive non-local Euclidean medians sparse unmixing (ANLEMSU) method is proposed. ANLEMSU improves NLSU by replacing the NLM total variation spatial consideration by the non-local Euclidean medians filtering approach, and it utilizes a joint maximum a posteriori (JMAP) (Mohammad-Djafari, 1996; Hsiao et al., 2002) strategy to acquire the relationships between the regularization parameters and the estimated fractional abundances. The contributions of this paper are twofold: (1) An efficient spatial consideration term that can cope with high-level noise, non-local Euclidean medians, is applied. This improves the traditional NLM method by reversing the top 50% of the weight values and abandoning the rest after sorting them in decreasing order and then computing the Euclidean medians with iteratively reweighted least squares (IRLS). (2) The regularization parameters are obtained adaptively, based on a JMAP strategy. In the ANLEMSU model, there are several unknown values: the objective fractional abundances and the two regularization parameters. To jointly estimate the unknowns is an ideal strategy. The ANLEMSU approach utilizes the relationships of the prior probability distributions between the abundances and the regularization parameters deduced from JMAP, and it estimates the regularization parameters and the unknown fractional abundances by an alternating iterative process. The rest of this paper is organized as follows. In Section 2, NLSU is reviewed. The ANLEMSU method is then proposed in Section 3, based on the techniques of non-local Euclidean medians (NLEM) and JMAP. Section 4 presents the experiments and analyses with two simulated datasets and two real hyperspectral remote sensing

images. Section 5 discusses the sensitivity of ANLEMSU in relation to the different parameters. Finally, the conclusions are drawn in Section 6. 2. Non-local sparse unmixing 2.1. Sparse unmixing In hyperspectral remote sensing imagery, each mixed pixel is represented as a vector y with L dimensions, and L is the number of spectral bands. Due to the computational tractability and flexibility, the LMM has been widely used in many applications, and the unmixing models discussed in this paper are all based on the LMM, shown as (1).

y ¼ Ma þ n

ð1Þ

where M is a L  q matrix containing q spectral signatures (called endmembers), and a is a q-dimensional vector containing the corresponding abundance of the endmember in M. n denotes the noise and model error and is also an L  1 vector. Because the components of a represent the fractional abundances, they should satisfy the abundance non-negative constraint (ANC) and abundance sumto-one constraint (ASC) (Heinz and Chang, 2001), as follows:

ai P 0 ði ¼ 1; 2; . . . ; qÞ q X

ð2Þ

ai ¼ 1

ð3Þ

i¼1

Sparse unmixing adopts another direction, which circumvents the determination of endmembers by selecting the optimal sparse combination of endmembers from a potentially quite large standard spectral library known in advance. Based on the LMM and the known spectral library, the sparse unmixing model is written as:

y ¼ Ax þ n

ð4Þ

where A acts as the large spectral library with a size of L  m, and x represents the m-dimensional vector as the abundance corresponding to library A, and is sparse. Similarly, the two constraints for the abundances should also be considered in the sparse unmixing model as:

minkxk0 s:t: ky  Axk2 6 d; x

x P 0;

1T x ¼ 1

ð5Þ

where ||x||0 represents the number of non-zero components in vector x, and d P 0 is the noise or model error. x P 0 and 1Tx = 1 denote the ANC and ASC, respectively. However, since the l0 term in (5) is a typical NP-hard problem, it can be relaxed as the l1 norm to obtain the sparsest solution under certain conditions (Candès and Romberg, 2007). Therefore, (5) can be replaced as follows:

minkxk1 s:t: ky  Axk2 6 d; x

x P 0;

1T x ¼ 1

ð6Þ

Pm

where kxk1 ¼ i¼1 jxi j, and xi represents the ith abundance in vector x. To tackle this convex optimization problem, SUnSAL was proposed. However, classical sparse unmixing focuses on analyzing the hyperspectral data without considering the spatial correlations. 2.2. Non-local sparse unmixing With the ongoing research into sparse unmixing techniques, the spatial information is now treated as important prior knowledge that can be incorporated into the conventional sparse unmixing model. NLSU was proposed based on an NLM method, combining the non-local spatial information (Manjón et al., 2008; Protter et al., 2009; Qian and Ye, 2013) and making systematic use of all possible self-predictions of the abundance maps. Zhong et al.

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

(2014a) illustrated the improvement of the NLSU over the classical sparse unmixing algorithm. However, although the use of the non-local similar windows of the abundance maps makes NLSU robust, the spatial information derived from the estimated abundances (or the weights which are used to represent the spatial correlation obtained from the intermediate estimated fractional abundances) is not so reliable. In the case of traditional NLM, the weights associated with each pixel are computed as follows:

  1 wij ¼ exp  2 kPi  Pj k2 h

ð7Þ

where wij denotes the weights of pixels i and j, and Pi and Pj represent the similar windows centered at pixels i and j, both locating in the non-local searching window. ||Pi  Pj||2 represents the squared l2 norm of the Euclidean distance between Pi and Pj, and h is a smoothing parameter (Buades et al., 2005a,b). In Fig. 1, we assume S is one searching window, and Pi and Pj are similar windows centered at pixels i and j, respectively. The final values computed by NLM are obtained by the weighted l2 norm of the distances, P NLM[xi] = jeS(i) wijxj, where xi and xj are the abundance values of x, located at pixel i and j. As two neighboring pixels with similar abundances for the same endmember will have a small Euclidean distance, the weight will be closer to 1. However, due to the inaccurate abundances, the NLSU algorithm can be further improved by dealing with the unreliable abundances. Furthermore, the regularization parameters of NLSU are usually kept fixed throughout the optimization process, and their values can significantly affect the unmixing performance. To obtain the optimal values of these parameters, NLSU often needs to be run multiple times with different parameter settings, which is time-consuming.

3. Adaptive non-local Euclidean medians sparse unmixing for hyperspectral imagery To better address the spatial information in sparse unmixing, ANLEMSU is proposed for hyperspectral remote sensing imagery, as an improvement of the NLM spatial consideration for sparse unmixing. However, differing from the NLSU algorithm, which considers the non-local spatial information by means of non-local total variation, ANLEMSU obtains the spatial correlations by means of filtering. Meanwhile, the regularization parameters are obtained adaptively, according to the relationships derived by JMAP.

11

3.1. The non-local Euclidean medians method Differing from the NLM method, which deals with each pixel of the abundance maps by means of averaging all the pixels whose neighborhood resembles the neighborhood of the current one, the NLEM method uses only part of the weights, together with their corresponding pixels, and achieves the final current pixel’s abundance by the use of IRLS (Chaudhury and Singer, 2012). In the non-local methods, both NLM and NLEM depend on the weights to associate the non-local spatial information with the current pixel. However, the high intensity of the inaccurate abundances existing in the non-local spatial domain leads to a negative influence. The NLEM method computes the weights with the original wij = exp(||Pi  Pj||2/h2) for every pixel j e S(i) in the searching window (S(i)) as (7). However, the final abundance of pixel i is determined by minimizing the summation of wij||Pi  Pj||, which is called the Euclidean median (or the geometric median) (Huber and Ronchetti, 2009), as the median is more robust to outliers in the weight values than the mean (Chaudhury and Singer, 2012, 2013). It should be noted that the NLM method actually uses the Euclidean mean over all the similar windows P, as the weights wij act as the averaged coefficient, which is equivalent to performing the following regression as (8) (Chaudhury and Singer, 2013). The NLM filtering process is to compute the minimizer of the weighted squared l2 norm of distance ||Pi  Pj||:

NLM½xi  ¼ arg min Pi

X

wij kPi  Pj k2

ð8Þ

j2SðiÞ

However, in the NLEM approach, the objective function is to solve the Euclidean median problem as (9) (Huber and Ronchetti, 2009) and the minimizer of the weighted l1 norm of these distances:

NLEM½xi  ¼ arg min Pi

X

wij kPi  Pj k

ð9Þ

j2SðiÞ

In practice, the NLEM method first reorders the weights in decreasing order, and the top 50% are then adopted, together with their corresponding abundances, to average:

P

j2SðiÞ Xð0Þ ¼ P

wij xj

j2SðiÞ

ð10Þ

wij

The initial Euclidean median for the current pixel is then achieved from (10). Since the idea of Euclidean medians is to use l1 regression instead of the l2 regression in NLM, to obtain a more reliable spatial correlation, the IRLS algorithm is employed. Consequently, the weights can be corrected gradually, according to the small bias v between the estimated Euclidean medians and the original neighborhood abundances.

P

ðkÞ

j2SðiÞ

XðkÞ ¼ P

wij vj xj ðkÞ

j2SðiÞ

wij vj

ð11Þ

where

1

vjðkÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðkxðk1Þ  xj k2 þ eÞ

ð12Þ

e is a small constant to prevent unstable results when ||x(k1)  xj||2 is close to zero. In addition, e was set as 1e4 in all the experiments Fig. 1. Scheme of the weight calculation in the non-local means strategy.

in this paper. The median x should be updated until the stopping criterion is satisfied, as ||x(k)  x(k1)||2 < e. Otherwise, k = k + 1.

12

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

3.2. Non-local Euclidean medians sparse unmixing (NLEMSU) Based on the non-local Euclidean medians method, the nonlocal Euclidean medians sparse unmixing algorithm (also called NLEMSU) is proposed as follows:

1 min kAX  Yk2F þ ksp kXk1;1 þ knl JNLEM ðXÞ þ iRmþ ðXÞ þ if1g ðXÞ X 2

ð13Þ

Formula (13) functions as a constraint optimization problem written as a version of minimizing the Lagrangian function, where Y e RLn represents the observed hyperspectral data matrix, and X e Rmn denotes the abundance matrix. A e RLm is the spectral library known in advance. In the NLEMSU model, kAX  Yk2F ¼ trððAX  YÞT ðAX  YÞÞ is the data fitting term, and P kXk1;1 ¼ nj¼1 kxj k1 acts as the sparse constraint, in which xj denotes the jth column of X. The spatial smoothness term written as JNLEM ðXÞ ¼ kNLEMðXÞ  Xk2F is used to incorporate the spatial relationships by means of the NLEM filtering approach in the abundance maps. The final terms iRqþ ðXÞ and i{1}(X) are the physical constraints considering the ground truth, which can be computed as:

iRmþ ðXÞ ¼ if1g ðXÞ ¼





0 minðxj Þ P 0; j ¼ 1; 2; . . . ; n þ1 minðxj Þ < 0; j ¼ 1; 2; . . . ; n P xj ¼ 1; P xj – 1; þ1

0

j ¼ 1; 2; . . . ; n

As the NLEMSU model (13) shows, there are two regularization parameters, ksp and knl , which need to be chosen properly. To realize the parameter adjustment adaptively, the NLEMSU problem is reconsidered by utilizing the JMAP technique, and the relationships between the regularization parameters and the estimated abundances can be obtained according to Bayes’ rule. In linear spectral unmixing problems, the aim is to find the typical endmembers’ signatures in the original image (Y) and their corresponding abundances (X) in each mixed pixel. However, in the case of the sparse unmixing approaches, the standard spectral library (A) can be obtained in advance. The problem of non-local spatial sparse unmixing can therefore be described as: Given Y and A, estimate X together with the parameters k = [ksp, knl]. As a prior, the abundance matrix X satisfies the sparse nature and follows a similar spatial distribution to the original observed image Y. Then, under the Bayesian framework, and following the main idea of JMAP (Mohammad-Djafari, 1996), the unmixing problem can be reconsidered as estimating the k parameters at the same level as estimating the abundance X. The original problem can be rewritten as:



b ksp ; knl ¼ arg max X;

ðX;ksp ;knl Þ

 pðX; ksp ; knl jYÞ

ð14Þ

Using Bayes’ rule, (14) can be expressed as:



b ksp ; knl Þ ¼ arg max pðYjX; ksp ; knl ÞpðXjksp ; knl Þpðksp ; knl Þ ð X; pðYÞ

ð15Þ

Because the estimated abundance map X and the regularization parameters k are independent of p(Y), (15) is equivalent to (16):

 b ksp ; knl Þ ¼ arg max pðYjX; ksp ; knl ÞpðXjksp ; knl Þpðksp ; knl Þ ð X;

!

ð16Þ

where p(Y|X, ksp, knl) is equivalent to p(Y|X) and is the likelihood distribution of the abundance map. Assuming that Y is contaminated with zero-mean white Gaussian noise, the probability distribution of p(Y|X) is as follows:

ð17Þ

where rn denotes the standard deviation of the Gaussian noise. For the abundance X, the prior model for the NLEM sparse unmixing is a combination of the Laplacian probability model for the sparse nature and a Gaussian distribution model for the spatial prior. In the ANLEMSU algorithm, the spatial prior information, addressed as the non-local Euclidean medians filtering, which accounts for the residual of the latest updated abundance X and the filtered one (NLEM(X)), can be viewed as obeying a general Gaussian probability distribution, as shown in (18):

! pffiffiffi 2  kXk1;1 kNLEMðXÞ  Xk2F 1 pðXjrX ; rres Þ ¼ exp   C1 rX 2r2res

ð18Þ

where C1 is a normalizing factor corresponding to the summation of the exponential function behind it. rX and rres denote the standard deviations of X and the residual image between the filtered abundance and the previous one (NLEM(X)  X), respectively. To some extent, estimating the regularization parameters is equivalent to estimating the standard deviations, as a result of their relationships ksp  (rX) and knl  (rres(NLEM(X)X)). The above prior models for X can then be rewritten as (19):

! pffiffiffi 2  kXk1;1 kNLEMðXÞ  Xk2F 1 pðXjksp ; knl Þ ¼ exp   C1 rX 2r2res

j ¼ 1; 2; . . . ; n

3.3. Adaptive non-local Euclidean medians sparse unmixing (ANLEMSU)



1 kY  AXk2F pðYjXÞ ¼ pffiffiffiffiffiffiffi exp  2r2n 2prn

ð19Þ

Furthermore, it should be noted that the prior for p(ksp, knl) is chosen to be a uniform distribution, as the values of ksp and knl are constants and range in a certain domain. To avoid under-flowing of the floating point numbers when computing the optimization problem, and to simplify the formula derived from the original MAP framework, the log function is brought to bear on (16) by plugging each proper prior introduced in (17) and (19). Then (16) yields: 

  b ksp ; knl ¼ arg max log PðYjXÞ þ log PðXjksp ; knl Þ þ logðksp ; knl Þ X;   8 9 <  1 2  kY  AXk2F þ log pffiffiffiffi1 = 2 rn 2prn   ¼ arg max pffiffi :  2  kXk  12  kNLEMðXÞ  Xk2 þ log 1 þ log C 2 ; 1;1 F C1 rX 2rres ( ) pffiffiffi 1 1 2 ¼ arg min  kY  AXk2F þ  kXk1;1 þ 2  kNLðXÞ  Xk2F 2 rX 2rn 2rres

ð20Þ where C2 refers to a normalized constant corresponding to the uniform prior of p(ksp, knl). The relationships between the unknown abundance and the regularization parameters can then be obtained from the above formulas as: pffiffi 8 2 < ^ksp ¼ r 2þren X 1

: ^k ¼ r2n nl 2r2res þe2

ð21Þ

where e1 and e2 are small constants used to prevent the denominators being close to zero, which would lead to unstable results. According to the relationships shown as (21), the values of the regularization parameters can be chosen when the abundance has been estimated. In addition, the method of estimating the above standard deviations should be mentioned here. The approach to estimating rn, which is called the average-filter method (Olsen, 1993), is simply to subtract the average filtered image from the noisy image, and then compute the noise at each pixel. To avoid some small-scale image structures from contributing to the noise measure, a small 3  3 local support window is chosen. The other standard deviation, rres, is obtained following the same method

13

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24 Table 1 SRE performance comparison of the estimated abundance maps. The bold values denote the higher precisions compared with others. Data

NNLS

SUnSAL

SUnSAL-TV

NLSU

NLEMSU

ANLEMSU

S-1

3.9294

4.0053 (ksps = 0.1)

13.644 (ksps = 1e5; ktv = 0.1)

14.3306 (ksps = 1e4; kspt = 300; mu = 15; h = 0.2)

16.1283 (ksps = 1e4; kspt = 300; mu = 14; h = 0.3)

16.1161 (mu = 8; h = 0.2)

S-2

4.2969

4.5724 (ksps = 0.1)

8.2779 (ksps = 0.03; ktv = 0.05)

9.9863 (ksps = 0.1; kspt = 300; mu = 0.9; h = 0.09)

11.0829 (ksps = 5e4; kspt = 10; mu = 2; h = 0.5)

10.0369 (mu = 8; h = 0.2)

used for estimating rn, after obtaining the residual image (NLEM(X)  X).

The process of how to compute the ANLEMSU model is introduced below in detail. Given the objective function (13), together with the adaptive update principles of the regularization parameters (21), the split augmented Lagrangian method of multipliers (ALM) and the alternating direction method of multipliers (ADMM) strategies (Eckstein and Bertsekas, 1992) are adopted. After inducing the certain intermediate variables for convenient representation, as shown as (22), the expanding of the objective function is developed as (23). 8 V1 > > > > < V2 s:t: V3 > > > V4 > : V5

1 kV1  Yk2F þ ksp kV2 k1;1 2

þ knl kV3 ðI  WÞk2F þ iRþ ðV4 Þ þ if1g ðV5 Þ þ

3.4. Scheme of computing the ANLEMSU model

1 min kV1 X;V1 ;V2 ;V3 ;V4 ;V5 2

LðX; V1 ; V2 ; V3 ; V4 ; V5 ; D1 ; D2 ; D3 ; D4 ; D5 Þ ¼

 Yk2F þ ksp kV2 k1;1 þ knl kV3 ðI  WÞk2F þ iRþ ðV4 Þ þ if1g ðV5 Þ

¼ AX ¼X ¼X ¼X ¼X ð22Þ

þ þ

l 2

l 2

kX  V2  D2 k2F þ

l 2

l 2

kX  V3  D3 k2F þ

kX  V5  D5 k2F

kAX  V1  D1 k2F

l 2

kX  V4  D4 k2F ð23Þ

where I is the identity matrix, and V3(I  W) represents the spatial consideration of the NLEM, which functions as V3  V3W by means of filtering. W is the weight matrix updated constantly by means of NLEM as soon as the latest estimated abundance matrix is acquired. l kAX  V1  D1 k2F and the other terms multiplied by the multiplier l2 2 (l > 0) are regarded as augmented terms related to the terms in the main function, such as kV1  Yk2F . The sequences D1, D2, D3, D4, D5 are also parts of the ADMM strategy satisfying: ( kþ1 kþ1 k kþ1 D1  AX þ V1 D1 . The main steps of comDkþ1 Dki  Xkþ1 þ Vkþ1 ði ¼ 2; 3; 4; 5Þ i i puting NLEMSU are listed in Algorithm 1.

Fig. 2. True fractional abundances of simulated dataset 1.

14

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

Algorithm 1: 1) Initialization: ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ set k=0, Max_Iter, kð0Þ sp ; knl , and l, and initially estimate X ; V 1 ; . . . ; V5 ; D1 ; . . . , and D5 , where 1

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

X0 ðAT A þ 4IÞ ðAT YÞ; V1 ¼ AXð0Þ , V2 ¼ Xð0Þ , V3 ¼ Xð0Þ ; V4 ¼ Xð0Þ , V5 ¼ Xð0Þ , D1 ¼ D2 ¼ D3 ¼ D4 ¼ D5 ¼ zerosðsizeðXÞÞ, and Y is the observed data. 2) Repeat: (2.1) The update of the regularization parameters:     ðkÞ kðkÞ f XðkÞ ; Y and knl f XðkÞ ; Y using (21); sp (2.2) The update of the abundance matrix x:   1  Xðkþ1Þ AT A þ 4I AT u1 þ u2 þ u3 þ u4 þ u5 where

u1 ¼ Vk1 þ Dk1 ; u2 ¼ Vk2 þ Dk2 ; u3 ¼ Vk3 þ Dk3 ; u4 ¼ Vk4 þ Dk4 ; u5 ¼ Vk5 þ Dk5 ; and I is the identity matrix. (2.3) The optimization of V1, V2, V3, V4, V5: ðkþ1Þ

(2.3.1)

V1

ðkþ1Þ

V1

ðkþ1Þ

(2.3.2)

V2

ðkþ1Þ

V2

 2  ðkÞ  arg minV1 12 kV1  Yk2F þ l2 AXðkÞ  V1  D1  F ; h  i ðkÞ ðkÞ 1  D1 1þl Y þ l AX  2 ðkÞ  l  ðkÞ arg minV2 kðkÞ  V 2  D2  sp kV2 k1 þ 2 X F   ; ðkÞ k soft u2 ; lsp

where soft(, s) functions as y ´ sign(y) max{|y|  s, 0}, as referred to in Zhang et al. (2010).  2  ðkþ1Þ ðkÞ ðkÞ  arg minV3 knl kV3 ðI  WÞk2F þ l2 XðkÞ  V3  D3  V3 F

ðkÞ 1  (2.3.3) (Zhang et al., 2010; Cai et al., 2009);  2knl ðkþ1Þ ðkÞ T ðkÞ V3 ðI  WÞðI  WÞ þ I X  D 3 l  2  ðkÞ  arg minV4 iRþ ðV4 Þ þ l2 XðkÞ  V4  D4  F ;   (2.3.4) ðkþ1Þ ðkÞ V4 ¼ max XðkÞ  D4 ; 0  2  ðkþ1Þ ðkÞ  arg minV5 if1g ðV5 Þ þ l2 XðkÞ  V5  D5  V5 F ;.   (2.3.5) ðkþ1Þ ðkÞ V5 ¼ XðkÞ  D5 þ ðH  JÞ h  i P where H = ones(n, 1) and J ¼ 1n 1  ni¼1 Xk  Dk5 . The  is denoted as the Kronecker product, and is defined as 2 3 a11 B    a1n B 6 .. .. .. 7 mn AB¼4 . , B e Rpq. . . 5 and A e R am1 B    amn B (2.4) Update the Lagrange multipliers: ðkþ1Þ

V4

ðkþ1Þ

(2.4.1) D1

ðkþ1Þ Di

ðkÞ

ðkþ1Þ

D1  AXðkþ1Þ þ V1 ðkÞ Di

ðkþ1Þ

;

ðkþ1Þ Vi

(2.4.2) X þ (i = 2, 3, 4, 5) 3) Update iteration: k k + 1; 4) Until stopping condition is satisfied and output the final abundance X.

In general, the stopping condition is set as a maximum number  2  2  2       of iterations (Max_Iter) or AXk  Vk1  þ Xk  Vk2  þ Xk  Vk3  þ F F F  2  2  k  k k k X  V4  þ X  V5  j < s, as the small s is a threshold such as F

F

s 6 1e5. In the process of ANLEMSU, the regularization parameters, ksp and knl, are initially chosen at random, and the abundance is estimated as X(0) (ATA + 4I)1(ATY). Consequently, the regularization parameters and the abundances are alternately updated with the given principle as (19), until the stopping criterion is satisfied. The computational complexity of the proposed ANLEMSU algorithm is as follows. We have assumed that the observed hyperspectral image Y is L  n, where L is the number of bands and n is the number of pixels. The spectral library A is L  m, and X e Rmn. Then, in Algorithm 1, the computational complexity for updating the regularization parameters (ksp and knl) is O(nL). In the step of computing X (step 2.2), the order of complexity is O(nL2), as proved in Iordache et al. (2012). For V1, V2, V4, V5 and D1, D2, D3, D4, D5, the

complexity is O(n). In the step of solving V3, the computational complexity of the update of the non-local weight matrix W is O(n2L) (Chaudhury and Singer, 2012, 2013), and the computational process of V3 is O(n3). Finally, taking all the computational complexity for updating X per iteration into account, the total computational complexity for the proposed algorithm is max(n2  (n, L)). It should be noted that n is more likely to be bigger than L, which leads to the overall complexity of ANLEMSU being O(n3), and the most costly steps in the proposed approach are the calculation of the spatial regularization term. 4. Experiments and analyses To evaluate the performance of the proposed algorithm, two simulated hyperspectral datasets and two real hyperspectral remote sensing images were used with the traditional spectral unmixing approach, non-negative least squares (NNLS), the classical sparse unmixing algorithm, SUnSAL, the spatial sparse unmixing techniques, SUnSAL-TV and NLSU, as well as NLEMSU and the

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

15

Fig. 3. True fractional abundances of simulated dataset 2.

proposed ANLEMSU. The accuracy assessments of all the experiments in this paper were made by the signal-to-reconstruction h i h i ^ k22 and SRE (dB) = 10 log10(error (SRE), SRE ¼ E kxk22 =E kx  x ^ SRE), where x is denoted as the true fractional abundance, and x is the estimated one (Iordache et al., 2011). In addition, the main parameters corresponding to each algorithm were set as shown in Table 1, following the SRE values. The sizes of the search window and the similar window in the non-local spatial sparse unmixing approaches of NLSU, NLEMSU, and ANLEMSU were all set as 5  5 and 3  3, respectively. 4.1. Experiments with simulated datasets In the first simulated experiment, the abundances of simulated dataset 1 (S-1) were generated by the LMM with 75  75 pixels for five endmembers, and the endmembers were randomly selected from the USGS spectral library, with 224 bands. In addition, the ASC and ANC were both imposed in the simulation process. There were both pure and mixed regions with different kinds and differ-

ent proportions of endmembers. The background was also highly mixed with all these five endmembers. Fig. 2 shows the true abundances of each of these endmembers, with their proportional scale on the right. Finally, a high level of Gaussian noise was added, with the SNR = 10 dB. The simulated data 2 (S-2), with 100  100 pixels and 221 bands per pixel, was created as a standardized simulated dataset for benchmarking the accuracy of the spectral unmixing techniques provided in the HyperMix tool (Jimenez et al., 2012). The dataset uses fractals and generates distinct spatial patterns across the image. As with S-1, the abundances of S-2 also satisfy the LMM, together with the two physical constraints. The endmembers for the simulation were randomly selected from a spectral library, which was built by the given signatures, together with the USGS library, after removing certain bands. Finally, zero-mean Gaussian noise was added with the SNR = 10 dB. The true fractional abundances of S-2 are shown in Fig. 3. The results obtained by these two simulated datasets are shown in Figs. 4 and 5, by the use of the NNLS, SUnSAL, SUnSAL-TV, NLSU, NLEMSU, and ANLEMSU algorithms, respectively. Assessments

16

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

Endmember 1

Algorithm

Endmember 2

Endmember 5

1 10

0.9 0.8

20

10

0.9 0.8

20

0.7 30

NNLS

0.6 0.5

40

0.4

50

0.3

60

0.2

70

0.1 10

20

30

40

50

60

70

30

0.6 0.5

40

0.4

50

0.3

60

0.2

70

0.1

0

0.9 0.8

20

10

20

30

40

50

60

SUnSAL

0.6 0.5

40

0.4

50

0.3

10

0.9 0.8

20 30

0.6 0.5

40

0.4

50

0.3 0.2

70

0.1

70

0.1

40

50

60

10

70

20

30

40

50

60

70

SUnSAL -TV

0.8

10

20

0.7

20

0.8 0.7 0.6

30

0.5 40

0.4

50

0.3 0.2

60

0 10

20

30

40

50

60

0.4

50

0.3 0.2

60

0.6

20

0.5

30

NLSU

0.4

40

0.3

50

0.2

60

0.1

70 10

20

30

40

50

60

70

70

0 10

0.7

0.7 0.6 20 0.5 30

20

30

40

50

60

0.3 50

20

0.6

30

0.5

20

30

40

50

60

70

10 20 30

ANLEMSU

0.4

40

0.3

50 0.2 60

20

30

40

50

60

70

0

0.4 0.3

60

0.2 0.1

70 10

20

30

40

50

60

70

0

0.9

10

0.8

20

0.7

30

0.6 0.5

40

0.4 50

0.3 0.2

60

0.1

70

0 20

30

40

50

60

70

0.9

10

0.8 20

0.7

30

0.6 0.5 0.4

0.2

60

0.1

70 10

20

30

40

50

60

70

0.2 0.1

70 10

0.9

0.7

20

0.3

60

0

0.8

10

20

30

40

50

60

70

10

0.8 0.7

30

0.6

0.5 0.4

0.4

0.3

50

60

0.2

60

10

20

30

40

50

60

70

10

0.7

20

0.3 0.2 0.1

70

0

0.8

0.5

40

50

0.1

10

20

30

40

50

60

70

0.5

0.8

20

0.7

0.4

40

50

0.3

50

60

0.2

60

10

20

30

40

50

60

70

0.6

30

40

0.1

0 0.9

10

0.6 30

0

20

0.6

30

70 10

0.7

50

50

0.1 70

0.8

0.5

0.3

0.7

0.5

0.9

0.6

0.4

0

0.6

0

40

50

70 10

70

30

40

0.1 70

60

20

40

0.2 60

50

0.7

40

40

40

1

10

0.8

0.4

NLEMSU

30

10

70

10

0

10

0.1 20

0.1

70

10

0.2

70

0.5 40

0.1 70

0.3

60

0 0.9

10

0.6

0.4

50

1

0.9

30

0.5

0.7

60

30

0.6

40

1

0.2

20

0.7

30

10

60

10

0.8

20

70

0.7 30

0.9

10

0.7

1

10

1

1

0.5 0.4 0.3 0.2 0.1

70

0

10

20

30

40

50

60

70

0

Fig. 4. Estimated abundances of endmembers 1, 2, and 5 for S-1.

were made from both the qualitative and quantitative aspects, and the parameter settings are listed along with the SRE values in Table 1.

It can be seen from Fig. 4 that the results obtained by NNLS and SUnSAL are both full of outliers, and a visual comparison of the results confirms the different effects of the different algorithms.

17

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

Endmember 3

Algorithm

NNLS

Endmember 5 1

20

0.9 20

30

0.8 30

1

30

40

0.7 40

0.8

40

0.6

50

0.5

60

0.4 0.3

80

0.2

100

20

40

60

80

1

0.2

40

40 0.6

50 60

0.4 0.2

20

40

60

80

100

1

20

40

60

80

100

0.8

50

0.6

60

0.4

70 80

10

0

100

1

10

20

20

30

0.8 30

40

40

0.8 0.7 0.6 0.5

60

0.4

70

0.3

80

0.2

90

0.1

100

20

40

60

80

100

0.2

20

40

60

80

100

0 1 0.9 0.8 0.7

0 1

10

0.9

20

0.8

30

90

90

0 1.2

10

0.8 30

80

0.7

40

0.6

50

0.5

60

0.4

70

0.3

80

0.2

90 100

0.1 20

40

60

80

100

0 1

10

0.9

20

0.8

30

0.7

0.6

40

0.6

50

0.5

50

0.5

60

0.4 60

0.4

60

0.4

70

70

0.3

70

0.3

80

0.2

0.6

50

80 90 100

0.2 80

0.2

90

0.1

0 20

40

60

80

100

100

1

0 20

40

60

80

100

0.1

90 100

1

0 20

40

60

80

100

1

10

0.9 10

0.9

10

0.9

20

0.8 20

0.8

20

0.8

30

0.7 30

0.7

30

0.7

40

0.6 40

0.6

40

0.6

50

0.5 50

0.5

50

0.5

60

0.4 60

0.4

60

0.4

70

0.3 70

0.3

70

0.3

80

0.2 80

0.2

80

0.2

90

0.1 90

0.1

90

100

NLEMSU

0.4

80

30

70

NLSU

0.6

70

20

0.9

50

60

20

1

10 20

50

100

100

10

100

1.2

0.1 90

90

SUnSAL -TV

10

10

70

SUnSAL

Endmember 7

1.1

20

40

60

80

100

100

0 1

20

40

60

80

100

0 1

100

0.1 20

40

60

80

100

0 1

10

0.9

10

0.9

10

20

0.8 20

0.8

20

30

0.7 30

0.7

30

40

0.6 40

0.6

40

50

0.5 50

0.5

50

0.5

60

0.4 60

0.4

60

0.4

0.3 70

0.3

70

0.3

0.2 80

0.2

80

0.2

0.1 90

0.1

90

70 80 90 100

20

40

60

80

100

0

100

0 20

40

60

80

100

100

0.9 0.8 0.7 0.6

0.1 20

40

60

80

100

0

1

ANLEMSU

10

0.9

20

0.8

30

0.7

0.6

40

0.6

0.5 50

0.5

50

0.5

60

0.4 60

0.4

60

0.4

70

0.3 70

0.3

70

0.3

80

0.2 80

0.2

80

0.2

90

0.1 90

0.1

90

10

0.9 10

0.9

20

0.8 20

0.8

30

0.7 30

0.7

40

0.6 40

50

100

20

40

60

80

100

0

100

20

40

60

80

100

0

100

0.1 20

40

60

80

100

0

Fig. 5. Estimated abundances of endmembers 3, 5, and 7 for S-2.

Compared with SUnSAL-TV and NLSU, NLEMSU can better suppress the wrong unmixing abundance values when the SNR = 10 dB. The inaccurate abundance values in SUnSAL-TV and NLSU have been almost corrected by NLEMSU. ANLEMSU obtains similar results to the NLEMSU algorithm. For example, the corners of the squares are difficult to unmix, due to the complex neighborhoods, and inaccurate abundances often appear in the traditional unmixing methods because of the low quality of the imagery, and the

abundance maps often reflect outliers as abnormally bright pixels. However, the fractional abundances achieved by NLEMSU and ANLEMSU can effectively remove the outliers, and therefore outperform NLSU. Fig. 5 shows a visual comparison of the fractional abundance maps estimated from NNLS, SUnSAL, SUnSAL-TV, NLSU, NLEMSU, and ANLEMSU. Compared with the results of NNLS and SUnSAL, which are full of outliers, SUnSAL-TV, NLSU, NLEMSU, and

18

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

5 10 15 20 25 30 35 40 45

5

(a) R-1

10

15

20

25

30

35

40

45

50

(b) Spectral library for R-1

(c) Reference data

(d) The approximate true abundance maps Fig. 6. The R-1 dataset.

250

Gray values

200

150

100

50

0 0

5

10

15

20

25

30

35

40

45

50

The number of bands

(a) R-2 imagery

(b) Spectral signatures

(c) Reference data

1 2 5 6 (d) Classification result

3 7

4 8

(e) The approximate true abundance maps Fig. 7. The R-2 dataset.

ANLEMSU perform well when the SNR = 10 dB. It can be seen that the abundance maps obtained from SUnSAL-TV and NLSU exhibit obvious blurred regions due to the existence of noise, while the results obtained by the NLEMSU and ANLEMSU algorithms are significantly closer to the true abundances, especially for the parts with high-density noise. The NLEMSU and ANLEMSU algorithms utilize the advantages of NLEM and effectively suppress the outliers by means of the Euclidean medians. Table 1 gives a quantitative assessment for NNLS, SUnSAL, SUnSAL-TV, NLSU, NLEMSU, and ANLEMSU, in terms of SRE for the two simulated datasets, S-1 and S-2. In the comparison of S-1, the sparse unmixing algorithms show an outstanding performance in the spectral unmixing. Compared with the SRE values of over 3.9 for NNLS and SUnSAL, SUnSAL-TV and NLSU reach as high as 14, which is a great improvement. In NLEMSU and ANLEMSU, the SRE values are even higher at more than 16. Due to the adaptive nature of ANLEMSU, which can save much energy in the parameter

choice, ANLEMSU is much more practical than NLEMSU. It can again be seen that ANLEMSU obtains similar results to NLEMSU. For S-2, NLEMSU improves the SRE from the 9.9863 obtained by NLSU to 11.0829, which is an improvement of 1.0966. The quantitative assessment is consistent with the visual comparison. In general, NLEMSU exhibits better spatial consideration and less outliers than NLSU when the SNR is low. 4.2. Experiments with real hyperspectral imagery The real hyperspectral datasets chosen for the testing were both acquired by a Nuance NIR imaging spectrometer (650–1100 nm and 10 nm spectral interval). One dataset was gathered in April 2012 and the other in November 2013. As these datasets were collected by the authors, the same scene on the same day (with a high spatial resolution for the ground truth) was also captured using a digital camera.

19

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

Algorithm

Component 1

1.6

10

1.4

15

1.2

20

NNLS

35

0.6

35

40

0.4

40

45

0.2 10

20

30

40

50

0 1.8 1.6

10

1.4

15

1.2

20

1

25

35

40

0.4

40

45

0.2

45

40

50

50

1.4

10

15

1.2

15

20

1

20

10

25

0.8 0.6

35

45 10

20

30

40

50

15 20 25

45

0

50

1

45 10

20

30

40

50

1.2

15

0.6

0.2

10

20

30

40

50

0

1

0.8

0.6

0

1

25

0.8

40 45 10

20

30

40

50

10

20

1

25

0.8

30

0.6

35

0.4

40 45

0.2

50

0

10

20

30

40

50

50

0

0.4

45 50

0.2 10

20

30

40

50

5

1.6

10

1.4 1.2

20

1

25

0.8 0.6

35

0.4

40 45 50

0.2 10

20

30

40

50

0

10 15

0.8

0.6

25

0.4

35 40

0.2

45 50

0 10

20

30

40

50 1

5 1

0.8

0.6

0.9

10

0.8

15

0.7

20

0.6

25

0.5 0.4

30 0.4

0.2

0.3

35

0.2

40

0.1

45 10

20

30

40

50

0

50

0 10

20

30

40

50

1.2

1

5 1

0.8

0.6

0.4

40

50

1.2

40

35

0

15

30

30

45

1.4

20

25

0.2

1.6

5

10

5

20

0.4

0.2

10

20

35

0.4

40

50

0.6

0.6

40

30

35

15

30

0.8

35

20

30

1.2

15

1

5

25

45

1.4

1.2

30 0.4

20

0.2

1.6

5 10

1.4

25

15 0.8

5 10

0.4

1.6

20

1.2

1.4

0.6

40

ANLEMSU

0.2

0.8

30

50

0 1.2

35 40

1

35

NLEMSU

50

25

0.4

1.6

5 10

50

40

30

30

40

NLSU

30

5

5

50

20

25

0.6

30

10

20

35

20

0.2

15

30

1.8

10

30 0.4

5

0.8

10

0.6

10

30

50

0.8

45 50

2

5

15

20

30

30

1

15

0.8

5

SUnSAL -TV

10

25

25

Component 3 1.2

5

1

50

SUnSAL

Component 2

1.8

5

0.2

0.9

10

0.8

15

0.7

20

0.6

25

0.5

30

0.4

35

0.3 0.2

40

0.1

45 10

20

30

40

50

0

1.2

5 10

1

15 0.8

20 25

0.6

30 0.4

35 40

0.2

45

50

0 10

20

30

40

50

5

0.9

10

0.8

15

0.7

20

0.6

25

0.5

30

0.4

35

0.3 0.2

40

0.1

45

50

0

10

20

30

40

50

0

50 10

20

30

40

50

Fig. 8. Estimated abundances of endmembers 1, 2, and 3 for R-1.

The first real hyperspectral dataset (R-1) is a 50  50 subset, with 46 bands in each pixel. The spectral library used for this data was collected from other Nuance datasets obtained by the same spectrometer during the same time period, and contains 52 pure materials (Fig. 6(b)). The ground truth for this data was obtained following the methodology of Xu et al. (2013), and the reference data are also shown in Fig. 6(c). The approximate true abundance maps are shown in Fig. 6(d).

Unlike R-1, the R-2 dataset contains slightly more complicated land cover, with 128  64 pixels and 46 bands. This dataset consists of more than eight kinds of land-cover types, including wood (1), fresh grass (2), dead leaves (3), gravel (4), background (5), shadow (6), ceramic-bowl (7), and acrylic shoelaces (8), as shown in Fig. 7(a), and the spectral signatures of these eight land-cover types are shown in Fig. 7(b). To allow a quantitative assessment, we also took the same scene by the use of a digital camera with

20

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

Endmember 1

Algorithm

Endmember 5

Endmember 8

NNLS 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.5

1

1.5

2

2.5

3

3.5

4

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.2

0.4

0.6

0.8

1.2

1

SUnSAL 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

1.1

0.6

0.7

SUnSAL -TV -0.2

0

0.2

0.4

0.6

0

0.8

0.5

1

1.5

2

-0.2 -0.1

0

0.1

0.2

0.3

0.4

0.5

NLSU -0.6

-0.4

-0.2

0

0.2

0.6

0.4

0.8

1

0

1

0

0.5

1

1.5

2

0

0.2

0.4

0.6

0

0.2

0.4

0.6

1

0.8

NLEMSU -0.4

-0.2

0

0.2

-0.4

-0.2

0

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

0.8

1

ANLEMSU 0.4

0.6

0.8

1

0

0.5

1

1.5

0

0.2

0.4

0.6

0.8

1

Fig. 9. Estimated abundances of endmembers 1, 5, and 8 for R-2.

384  192 pixels (for its high spatial resolution) as the reference data, as shown in Fig. 7(c). Geometrical calibration, classification, and down-sampling were then undertaken on the reference data. The results of the classification are shown in Fig. 7(d), which includes the eight major land-cover types listed above, along with the background, and the approximate true abundance images were obtained after down-sampling of the classification results, as shown in Fig. 7(e). The classification process was undertaken in the ENVI 5.0 platform by the use of an objected-oriented method. The performances of the proposed NLEMSU and ANLEMSU, together with NNLS, SUnSAL, SUnSAL-TV, and NLSU, were tested using the R-1 and R-2 datasets. Figs. 8 and 9 show the abundance maps estimated for these two datasets, respectively. Table 2 lists the quantitative assessments for the two datasets. From Figs. 8 and 9, it can be seen that the abundance maps obtained by NLSU, NLEMSU, and ANLEMSU look alike visually. Compared with the NNLS method, SUnSAL obtains sparse results

and illustrates clear borders between the different blocks, especially for R-1. Unlike the plain distribution of the abundance maps acquired by NLSU, the abundance maps obtained by NLEMSU and ANLEMSU are more smooth and homogeneous. As can be seen in the details in component 2 of R-1, some outliers are suppressed, and more values with close abundances, which are probably inliers, are preserved. In R-2, the abundance map of the background, denoted as endmember 5, is quite close to the abundances obtained by the NLEMSU and ANLEMSU approaches. As a result of the medians avoiding the influence of outliers, the final abundances can achieve much more smooth characterizations of the material. Table 2 displays a quantitative comparison of NNLS, SUnSAL, SUnSAL-TV, NLSU, NLEMSU, and ANLEMSU for the R-2 dataset. Compared with NNLS, SUnSAL, and SUnSAL-TV, the SRE values of NLSU, NLEMSU, and ANLEMSU are a little higher, being over 6.00 for R-1 and close to 3.5 for R-2. The reason for this mainly lies in the different consideration of the spatial information. The

21

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24 Table 2 SRE performance comparison for the nuance estimated abundance maps. The bold values denote the higher precisions compared with others. Data

NNLS

SUnSAL

SUnSAL-TV

NLSU

NLEMSU

ANLEMSU

R-1

4.377

4.9281 (ksps = 500)

5.0392 (ksps = 200; ktv = 500)

6.0017 (ksps = 1e3; kspt = 300; mu = 140; h = 0.5)

6.0298 (ksps = 1e1; kspt = 300; mu = 150; h = 0.1)

6.0294 (mu = 150; h = 0.1)

R-2

1.3488

1.4692 (ksps = 240)

2.2645 (ksps = 2e5; ktv = 300)

3.0663 (ksps = 5e2; kspt = 5; mu = 150; h = 5)

3.2747 (ksps = 5e3; kspt = 300; mu = 160; h = 17)

3.4015 (mu = 190; h = 5)

8.280

8

7.457

7

6.635

10.00

10

8.856 7.713

8

5.813

6.569

4.990 4.167

5

3.345

4

4.281 3.138

4

2.523

3

5.425

6

SRE

SRE

6

1.700

1.994 0.8500

2

2 05 1E-

3 1E-

1E4

10 30 0 0

10 30 0 0

1E3

10

1

ps

10

0.1 0.0 1

0. 1

pt

0. 1

λs

ps

λs

1E3

1E5

3 1E-

0. 01

0.1

0.0 1

10

(a) SRE in relation to λsps and λspt for SUnSAL-TV

300 100

1

pt

4 1E-

λs

3 10000

1

10

1

λs

4 1E-

0. 01

-5 1E1

1E4

1E5

(b) SRE in relation to λsps and λspt for NLSU 11.10 9.669

10

8.238 6.806

8

SRE

5.375 6

3.944 2.513

4

1.081 -0.3500

2 0 -5 1E

10

300 100

0.1 0.0 1

0. 1

10 30 0 0

10

1

ps

λs

1E3

λs

-3 1E

pt

1

0. 01

-4 1E

1E4

1E -5

(c) SRE in relation to λsps and λspt for NLEMSU Fig. 10. SRE in relation to ksps and kspt for the S-2 dataset.

non-local spatial sparse unmixing approaches seem to be effective in improving the unmixing accuracy. In particular, the NLEM method for spatial correlation demonstrates an advantage over NLM for sparse unmixing. In our future research, more real hyperspectral datasets will be utilized to further test the effectiveness of the NLEMSU method. 5. Sensitivity analysis 5.1. Impact of the regularization parameters in the traditional algorithms For the spatial sparse unmixing algorithms, SUnSAL-TV, NLSU, and NLEMSU, one of the most important tasks is to determine the values of the regularization parameters ksps (the regularization parameter connected with the sparsity regularization term) and kspt (the parameter connected with the spatial constraint term). To analyze the impact of the regularization parameters in these

algorithms, and to illustrate the best parameter values, we plotted the curves of the SRE values with a series of different combinations of ksps and kspt, as shown in Fig. 10. The values of the regularization pairs for SUnSAL-TV, NLSU, and NLEMSU were set as follows: ksps and kspt were both set as {1e5 5e5 1e4 5e4 1e3 5e3 1e2 3e2 5e2 1e1 1 10 100 300}. Simulated dataset S-2 with the SNR = 10 dB was used as an example to better illustrate the impact of the regularization parameters. The relationship between the SREs and the two regularization parameters is shown in Fig. 10, which represents the searching process for the optimal combinations of regularization parameters. It can be seen that a proper value of each prior term can lead to a better unmixing accuracy. The curves of Fig. 10 show part of the results from the spatial sparse unmixing computation. It is actually quite time-consuming to obtain the best regularization parameters, which means that the corresponding adaptive methods are necessary. Given the trend of these curves, it can be seen that

22

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

9.990

10.0

9.646

10.51

9.303 8.959 8.615

9.0

8.271 7.928

8.5

9.980 9.450

10

8.920

SRE

9.5

8.390

9

7.860

7.584

7

0.0 9

0.0 8

0.2

1

μ

4

0.0 1

5

0. 9

0.4

3

0.0 2

0. 8

μ

0. 6

0.6

h

0. 4

0.9

0.8

0.5

0.0 6

0.0 4

6.800

(a) SRE in relation to mu and h for NLSU

0.1

h

7.5 0.1

7.330 8

7.240

8.0

2

SRE

11.04 11

0.2

(b) SRE in relation to mu and h for NLEMSU 10.05

10

8.888 7.725

8

6.563

SRE

5.400

6

4.238 3.075

4

1.913 0.7500

2 0.8

2

0.6

6

0.4

h

4

0

μ 10

8

0.2

(c) SRE in relation to mu and h for A-NLEMSU Fig. 11. SREs in relation to mu and h for S-2.

the sparse regularization parameter ksps should be small, in general, compared with the spatial parameter kspt, e.g., kspt = 10, but ksps = 5e4 for NLEMSU. 5.2. Sensitivity to the inner parameters of the non-local methods In the non-local spatial sparse unmixing approaches, apart from the two regularization parameters, ksps and kspt, there are also two inner parameters in the non-local method, mu and h. mu, which

acts as a Lagrangian multiplier, decides the degree of split in the optimization procedure; and h is the decay coefficient, controlling the filtering degree in the non-local methods. The values of mu and h significantly influence the accuracy of the unmixing. It is therefore essential to obtain the appropriate values in the non-local spatial sparse unmixing methods. Similarly, S-2 was used for the sensitivity analysis of mu and h. The optimal values of mu and h were searched for over the following ranges: mu = {0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9} or {1 2 3

Fig. 12. SRE performance comparison for S-2 with different SNRs.

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

4 5 6 7 8 9} and h={0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009} or {0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09}. The relationships between mu, h, and the SREs of S-2 are displayed in Fig. 11. As shown in Fig. 11, different combinations of mu and h lead to different SREs for the non-local spatial sparse unmixing methods. It can also be seen that the values of h are correlated with the imagery, due to the same noise or similar unmixing errors, e.g., h is close to 0.1 for all the non-local approaches in the S-2 dataset. However, the values of mu mostly depend on the spatial prior term itself, and it accounts for the intermediate spatial relationship during the iterative process. In summary, the non-local inner parameters do have certain effects on the final results.

23

Acknowledgements This work was supported by National Basic Research Program of China (973 Program) under Grant No. 2011CB707105, the National Natural Science Foundation of China under Grant No. 41371344, the 863 High Technology Program of the People’s Republic of China under Grant No. 2013AA12A301, A Foundation for the Author of National Excellent Doctoral Dissertation of the People’s Republic of China (FANEDD) under Grant No. 201052, and the Fundamental Research Funds for the Central Universities under Grant No. 2042014kf00231, and Program for Changjiang Scholars and Innovative Research Team in University (IRT1278). The authors would like to thank the editor, associate editor and anonymous reviewers for their helpful comments and suggestions.

5.3. Sensitivity of the sparse unmixing algorithms to the SNR References Different SNRs means different image qualities. To test the sensitivity of the sparse unmixing algorithms to the SNR, we added one more case of simulated dataset S-1 with a higher SNR, which was 20 dB. The higher the quality of the image is, the easier the unmixing process should be. The quantitative assessment for this added dataset, together with the aforementioned S-2 with the SNR = 10 dB, is shown in Fig. 11. In Fig. 12, the left (red1) group is the results of the S-2 dataset with the SNR = 10 dB, and the right (blue) group is the results with the SNR = 20 dB. It can be seen that NLEMSU and ANLEMSU show a greater superiority with a lower SNR, which means that the nonlocal Euclidean median sparse unmixing methods are more robust with regard to strong noise and the influence of outliers. As the quality of remote sensing imagery cannot be guaranteed, due to the effects of the complicated environment during the data acquisition and transmission, the better anti-outlier capability makes the proposed approach much more practical in reality.

6. Conclusion This paper proposes a new adaptive spatial sparse unmixing algorithm based on a non-local Euclidean medians filtering method (ANLEMSU). To better illustrate the effectiveness of ANLEMSU, the traditional spectral unmixing method, NNLS, together with the latest sparse and spatial sparse unmixing techniques, SUnSAL, SUnSAL-TV, and NLSU, were used to confirm the advantages of the novel algorithms, NLEMSU and the proposed ANLEMSU. Differing from non-local sparse unmixing (NLSU), which uses the non-local means (NLM) total variational method for the spatial information consideration, non-local Euclidean medians (NLEM) is incorporated into the sparse unmixing model by means of filtering. Due to the complete consideration of the non-local spatial information in NLEM, NLEMSU and ANLEMSU achieve better unmixing results with both simulated datasets and real hyperspectral images. In addition, ANLEMSU can adjust the regularization parameters adaptively, based on the JMAP estimation technique, which makes ANLEMSU more practical in real applications. The experimental results with both simulated data and real hyperspectral images confirm the superiority of NLEMSU and ANLEMSU, especially in coping with outliers in low-quality images, since the NLEM method abandons 50% of the less reliable weights, together with the corresponding similar windows, and computes the final abundance values by the IRLS algorithm. Our future research will focus on improving the computational efficiency when computing the non-local Euclidean medians weight matrix. 1 For interpretation of color in ‘Fig. 12’, the reader is referred to the web version of this article.

Bioucas-Dias, J.M., Plaza, A., Dobigeon, N., Parente, M., Du, Q., Gader, P., Chanussot, J., 2012. Hyperspectral unmixing overview: geometrical, statistical, and sparse regression-based approaches. IEEE J. Sel. Top. Appl. Earth Observ. Rem. Sens. (JSTARS) 5 (2), 354–379. Bioucas-Dias, J.M., Plaza, A., Camps-Valls, G., Scheunders, P., Nasrabadi, N.M., Chanussot, J., 2013. Hyperspectral remote sensing data analysis and future challenges. IEEE Geosci. Rem. Sens. Mag. 1 (2), 6–36. Buades, A., Coll, B., Morel, J.M., 2005a. A review of image denoising algorithms, with a new one. Multiscale Model. Simul. 4 (2), 490–530. Buades, A., Coll, B., Morel, J.M., 2005. A non-local algorithm for image denoising. In: IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR 2005), vol. 2, pp. 60–65. Cai, J.F., Osher, S., Shen, Z., 2009. Convergence of the linearized Bregman iteration for l1-norm minimization. Math. Comput. 78, 2127–2136. Candès, E., Romberg, J.J., 2007. Sparsity and incoherence in compressive sampling. IEEE Trans. Image Process. 23, 969–985. Chaudhury, K.N., Singer, A., 2012. Non-local Euclidean medians. IEEE Signal Process. Lett. 19 (11), 745–748. Chaudhury, K.N., Singer, A., 2013. Non-local patch regression: robust image denoising in patch space. In: IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 1345–1349. Demarchi, L., Chan, J.C.-W., Ma, J., Canters, F., 2012. Mapping impervious surfaces from superresolution enhanced CHRIS/Proba imagery using multiple endmember unmixing. ISPRS J. Photogram. Rem. Sens. 72, 99–112. Eckstein, J., Bertsekas, D., 1992. On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators. Math. Program. 55, 293–318. Heinz, D.C., Chang, C.-I., 2001. Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery. IEEE Trans. Geosci. Rem. Sens. 39 (3), 529–545. Hsiao, I.-T., Rangarajan, A., Gindi, G., 2002. Joint-MAP Bayesian tomographic reconstruction with a Gamma-mixture prior. IEEE Trans. Image Process. 11 (12), 1466–1477. Huber, P.J., Ronchetti, E.M., 2009. Robust statistics. In: Wiley Series in Probability and Statistics. Wiley, Hoboken, NJ. Iordache, M.D., 2011. A sparse regression approach to hyperspectral unmixing. Ph.D. Dissertation, School of Elect. Comput. Eng., IST. Iordache, M.D., Bioucas-Dias, J.M., Plaza, A., 2011. Sparse unmixing of hyperspectral data. IEEE Trans. Geosci. Rem. Sens. 49 (6), 2014–2039. Iordache, M.D., Bioucas-Dias, J.M., Plaza, A., 2012. Total variation spatial regularization for sparse hyperspectral unmixing. IEEE Trans. Geosci. Rem. Sens. 50 (11), 4484–4502. Iordache, M.D., Bioucas-Dias, J.M., Plaza, A., 2014a. Collaborative sparse regression for hyperspectral unmixing. IEEE Trans. Geosci. Rem. Sens. 52 (1), 341–354. Iordache, M.D., Bioucas-Dias, J.M., Plaza, A., Somers, B., 2014b. MUSIC-CSR: hyperspectral unmixing via multiple signal classification and collaborative sparse regression. IEEE Trans. Geosci. Rem. Sens. 52 (7), 4364–4382. Jimenez, L.I., Martin, G., Plaza, A., 2012. A new tool for evaluating spectral unmixing applications for remotely sensed hyperspectral image analysis. In: Proceeding of the 4th International Conference on Geographic Object-Based Image Analysis (GEOBIA), Rio de Janeiro, Brazil, pp. 1–5. Keshava, N., Mustard, J.F., 2002. Spectral unmixing. IEEE Signal Process. Mag. 19 (1), 44–57. Li, C., Fang, F., Zhou, A., Zhang, G., 2014. A novel blind spectral unmixing method based on error analysis of linear mixture model. IEEE Geosci. Rem. Sens. Lett. 11 (7), 1180–1184. Liu, J., Zhang, J., 2014. Spectral unmixing via compressive sensing. IEEE Trans. Geosci. Rem. Sens., 10.1109/TGRS.2014.2307573. Ma, W.K., Bioucas-Dias, J.M., Tsung-Han, C., Gillis, N., Gader, P., Plaza, A.J., Ambikapathi, A., Chong-Yung, C., 2014. A signal processing perspective on hyperspectral unmixing: Insights from remote sensing. IEEE Signal Process. Mag. 31 (1), 67–81.

24

R. Feng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 97 (2014) 9–24

Manjón, J.V., Carbonell-Caballero, J., Lull, J.J., García-Martí, G., Martí-Bonmatí, L., Robles, M., 2008. MRI denoising using non-local means. Med. Image Anal. 12 (4), 514–523. Mohammad-Djafari, A., 1996. Joint estimation of parameters and hyperparameters in a Bayesian approach of solving inverse problems. In: International Conference on Image Processing, vol. 2, pp. 473–476. Olsen, S.I., 1993. Noise variance estimation in images. In: Proceeding of the 8th Scandinavian Conference on Image Analysis (SCIA), Tromso, Norway, pp. 25–28. Plaza, A., Du, Q., Bioucas-Dias, J.M., Jia, X., Kruse, F., 2011. Foreword to the special issue on spectral unmixing of remotely sensed data. IEEE Trans. Geosci. Rem. Sens. 49 (11), 4103–4110. Protter, M., Elad, M., Takeda, H., Milanfar, P., 2009. Generalizing the nonlocalmeans to super-resolution reconstruction. IEEE Trans. Image Process. 18 (1), 36–51. Pu, H., Xia, W., Wang, B., Jiang, G., 2014. A fully constrained linear spectral unmixing algorithm based on distance geometry. IEEE Trans. Geosci. Rem. Sens. 52 (2), 1157–1176. Qian, Y., Ye, M., 2013. Hyperspectral imagery restoration using nonlocal spectralspatial structured sparse representation with noise estimation. IEEE J. Sel. Top. Appl. Earth Observ. Rem. Sens. (JSTARS) 6 (2), 499–515. Qu, Q., Nasrabadi, N.M., Tran, T.D., 2014. Abundance estimation for bilinear mixture models via joint sparse and low-rank representation. IEEE Trans. Geosci. Rem. Sens. 52 (7), 4404–4423. Tang, W., Shi, Z., Duren, Z., 2014. Sparse hyperspectral unmixing using an approximate L0 norm. Optik-Int. J. Light Electron Opt. 125 (1), 31–38.

Tits, L., Keersmaecker, W.D., Somers, B., Asner, G.P., Farifteh, J., Coppin, P., 2012. Hyperspectral shape-based unmixing to improve intra- and interclass variability for forest and agro-ecosystem monitoring. ISPRS J. Photogram. Rem. Sens. 74, 163–174. Tong, Q., Xue, Y., Zhang, L., 2014. Progress in hyperspectral remote sensing science and technology in China over the past three decades. IEEE J. Sel. Top. Appl. Earth Observ. Rem. Sens. (JSTARS) 7 (1), 70–91. Wu, Z., Wei, Z., Liu, J., 2011. A semi-supervised hyperspectral unmixing method based on enhanced sparsity. Proc. Eng. 15, 1864–1868. Xu, X., Zhong, Y., Zhang, L., Zhang, H., 2013. Sub-pixel mapping based on a MAP model with multiple shifted hyperspectral imagery. IEEE J. Sel. Top. Appl. Earth Observ. Rem. Sens. (JSTARS) 6 (2), 580–593. Zhang, X., Burger, M., Bresson, X., Osher, S., 2010. Bregmanized nonlocal regularization for deconvolution and sparse reconstruction. SIAM J. Imaging Sci. 3 (3), 253–276. Zhao, X., Wang, F., Huang, T., Ng, M.K., Plemmons, R.J., 2013. Deblurring and sparse unmixing for hyperspectral images. IEEE Trans. Geosci. Rem. Sens. 51 (7), 4045– 4058. Zhong, Y., Feng, R., Zhang, L., 2014a. Non-local sparse unmixing for hyperspectral remote sensing imagery. IEEE J. Sel. Top. Appl. Earth Observ. Rem. Sens. (JSTARS) 7 (6), 1889–1909. Zhong, Y., Zhao, L., Zhang, L., 2014b. An adaptive differential evolution endmember extraction algorithm for hyperspectral remote sensing imagery. IEEE Geosci. Rem. Sens. Lett. 11 (6), 1061–1065. Zhu, F., Wang, Y., Xiang, S., Fan, B., Pan, C., 2014. Structured sparse method for hyperspectral unmixing. ISPRS J. Photogram. Rem. Sens. 88, 101–118.