IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 51, NO. 7, JULY 2013 4045 Deblurring and Sparse Unmixing for Hyperspectral Images Xi-Le Zhao, Fan Wang, Ting-Zhu Huang, Michael K. Ng, and Robert J. Plemmons Abstract—The main aim of this paper is to study total vari- ation (TV) regularization in deblurring and sparse unmixing of hyperspectral images. In the model, we also incorporate blurring operators for dealing with blurring effects, particularly blurring operators for hyperspectral imaging whose point spread functions are generally system dependent and formed from axial optical aberrations in the acquisition system. An alternating direction method is developed to solve the resulting optimization problem efficiently. According to the structure of the TV regularization and sparse unmixing in the model, the convergence of the alternating direction method can be guaranteed. Experimental results are reported to demonstrate the effectiveness of the TV and spar- sity model and the efficiency of the proposed numerical scheme, and the method is compared to the recent Sparse Unmixing via variable Splitting Augmented Lagrangian and TV method by Iordache et al. Index Terms—Alternating direction methods, deblurring, hyperspectral imaging, linear spectral unmixing, total variation (TV). I. I NTRODUCTION H YPERSPECTRAL sensors collect multiple images of a scene using hundreds of contiguous spectral bands from ultraviolet to visible to infrared. There are a wide range of applications such as remote surveillance, military target discrimination, medicine, astrophysics, and environmen- tal monitoring, e.g., [1] and [2]. Many data analysis methods have been proposed and developed for hyperspectral images, for instance, classification, clustering, and spectral unmixing Manuscript received May 19, 2012; revised October 6, 2012; accepted November 1, 2012. Date of publication January 29, 2013; date of current version June 20, 2013. The work of M. K. Ng was supported by the Hong Kong Research Grant Council and Hong Kong Baptist University Faculty Research Grant. The work of X.-L. Zhao and T.-Z. Huang was supported in part by the National Natural Science Foundation of China under Grant 61170311, by the Chinese Universities Specialized Research Fund for the Doctoral Program under Grant 20110185110020, and by the Sichuan Province Science and Technology Research Project under Grants 2011JY0002 and 12ZC1802. The work of R. J. Plemmons was supported in part by the U.S. Air Force Office of Scientific Research, under Grant FA9550-11-1-0194, and in part by the U.S. National Geospatial-Intelligence Agency under Contract HM1582-10-C-0011. (Corresponding author: M. K. Ng.) X.-L. Zhao and T.-Z. Huang are with the School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu 610054, China (e-mail: xlzhao122003@163.com; tingzhuhuang@126.com). F. Wang is with the Department of Mathematics and Statistics, Lanzhou University, Lanzhou, China (e-mail: 09466029@hkbu.edu.hk). M. K. Ng is with the Centre for Mathematical Imaging and Vision and the Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong (e-mail: mng@math.hkbu.edu.hk). R. J. Plemmons is with the Department of Computer Science and the De- partment of Mathematics, Wake Forest University, Winston-Salem, NC 27106 USA (e-mail: plemmons@wfu.edu). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2012.2227764 [3]–[8]. In these methods, one is often interested in determining the underlying materials in each pixel. By assuming that the measured spectrum of each mixed pixel is a linear combination of spectral signatures (called endmembers), the linear spectral mixture model can be formulated as follows: Y = XA + N (1) where the matrix A ∈ R m×l (m is usually larger than l) repre- sents a spectral library containing m spectral signatures with l spectral bands, Y ∈ R n×l is the observed data matrix (each row contains the observed spectrum of a given pixel), X ∈ R n×m is the fractional abundance matrix of endmembers (each column contains the fractional abundances of a given endmember), and N ∈ R n×l is the matrix collecting the errors or noise affecting the measurements at each spectral band. In the following dis- cussion, we denote the (i, j )th entry of a matrix M by M i,j or (M ) i,j . Given the matrices Y and A, our purpose is to find the endmembers present in pixels, i.e., to determine the matrix X. In the literature, the abundance nonnegativity constraint X ≥ 0 (i.e., each entry is nonnegative) and the abundance sum-to-one constraint m  j=1 X i,j =1, i =1,...,n are generally imposed on the fractional abundances [8]. There are several traditional methods used for such unmix- ing, including nonnegatively constrained least squares method, pure-pixel-based algorithms (e.g., N-FINDR [9]), minimum- volume-based algorithms (e.g., Simplex Identification via Split Augmented Lagrangian [10]), and statistical algorithms [11] to determine the fractional abundances; see [8] for a detailed discussion. Recently, several new unmixing methods based on compressed sensing have been developed and studied, see, e.g., [6], [12], [17], [19], [20], [23], [24], and [26]. Gillis and Plemmons [17], [18] have proposed the use of nonnegative matrix underapproximation to extract features in a recursive way so that hyperspectral data can be com- pressively represented. They also used a sparsity constraint to extract materials in each pixel [18]. In [19], Guo et al. studied l 1 unmixing model and developed fast computational approaches based on Bregman iterations. In [12], Bioucas- Dias and Figueiredo studied alternating direction algorithms for constrained sparse regression for hyperspectral unmixing. In [6], Iordache et al. studied and compared several available and new linear sparse regression algorithms to solve the spectral 0196-2892/$31.00 © 2013 IEEE