自动化学报  2018, Vol. 44 Issue (1): 116-128   PDF    
一种基于协同稀疏和全变差的高光谱线性解混方法
陈允杰1, 葛魏东1, 孙乐2     
1. 南京信息工程大学数学与统计学院 南京 210044;
2. 南京信息工程大学计算机与软件学院 南京 210044
摘要: 稀疏分解是高光谱图像(Hyperspectral image,HSI)解混中的常用方法,为了克服传统稀疏解混方法只重视挖掘空间相关性而忽视稀疏性精确刻画的缺点,本文提出一种新的基于协同稀疏和全变差(Total variation,TV)相结合的高光谱空谱联合线性解混方法,从而进一步提高解混的精度.该方法基于已知光谱库的高光谱稀疏线性回归模型,利用TV正则项对高光谱邻域像元间的相关性进行约束;同时,协同稀疏性被用来刻画丰度系数的行稀疏性,从而表明协同稀疏先验对空谱联合解混精度的提高至关重要;最后采用交替方向乘子法求解模型.模拟高光谱数据实验结果定量地验证本文方法能够比现有同类方法获得更精确的解混结果,同时真实高光谱数据实验结果定性地验证了本文方法的有效性.
关键词: 高光谱图像     协同稀疏     TV正则项     线性光谱解混     交替方向乘子法    
A Novel Linear Hyperspectral Unmixing Method Based on Collaborative Sparsity and Total Variation
CHEN Yun-Jie1, GE Wei-Dong1, SUN Le2     
1. College of Math and Statistics, Nanjing University of Information Science and Technology, Nanjing 210044;
2. College of Computer and Software, Nanjing University of Information Science and Technology, Nanjing 210044
Manuscript received : May 23, 2016, accepted: December 10, 2016.
Foundation Item: Supported by National Natural Science Foundation of China (61672291, 61601236, 61471199, 61571230) and Natural Science Foundation of Jiangsu Province (BK20150923)
Author brief: CHEN Yun-Jie  Professor at the College of Math and Statistics, Nanjing University of Information Science and Technology. He received his Ph. D. degree from Nanjing University of Science and Technology in 2008. His research interest covers image processing, pattern recognition, and numerical analysis;
GE Wei-Dong  Master student at the College of Math and Statistics, Nanjing University of Information Science and Technology. His main research interest is image processing of hyperspectral unmixing
Corresponding author. SUN Le  Lecturer at the School of Computer and Software, Nanjing University of Information Science and Technology. He received his Ph. D. degree from Nanjing University of Science and Technology in 2014. His research interest covers image processing especially for hyperspectral imagery, including unmixing, classification, and restoration. Corresponding author of this paper
Recommended by Associate Editor XIN Jing-Min
Abstract: Sparse decomposition is one of the popular tools for hyperspectral unmixing. In order to overcome the shortcomings of traditional sparse unmixing methods which only pay attention to the spatial correlation and neglect depicting sparsity accurately, we propose a new spatial-spectrally linear hyperspectral unmixing method based on collaborative sparsity and total variation (TV) regularization to further improve the accuracy of unmixing. This method is based on hyperspectral sparse linear regression model with a spectral library given in advance, in which the total variation is utilized to impose a constraint on the correlation between neighboring pixels of hyperspectral image (HSI). Meanwhile, the collaborative sparsity is explored to depict the row-sparse characteristic of the fractional abundances, thus pointing out the fact that the collaborative sparsity prior plays an important role in further accuracy improvement of HSI spatial-spectral unmixing. At last, the proposed model is solved by the alternating direction method of multipliers. Experimental results on simulated hyperspectral data quantitatively validate that the our method outperforms those state-of-the-art algorithms, and the experimental results on real hyperspectral data qualitatively verify the effectiveness of the algorithm.
Key words: Hyperspectral image (HSI)     collaborative sparsity     total variation (TV)     linear spectral unmixing     alternating direction method of multipliers    

20世纪80年代以来, 高光谱遥感技术发展显著.由于其具有较高的光谱分辨率, 高光谱遥感已广泛应用于地质制图、植被调查、大气研究、环境监测、目标侦察、伪装识别等重要领域[1].

由于光谱成像仪的限制和地物的复杂多样性, 高光谱传感器获得的图像中往往包含多种物质(端元)的混合像元[2-3].为了提高后续应用的精度, 从混合像元中分解与提取各个地物光谱并求得其对应比例至关重要, 即将混合像元分解成不同的端元组合, 并求得端元在该像元中所占的比例, 即丰度系数.

线性光谱混合模型是高光谱解混常用的模型[4-5], 其假设每个像元可由端元线性表示.基于稀疏性的线性高光谱解混方法可分为基于稀疏性约束的非负矩阵分解的解混方法和基于已知光谱库的稀疏解混方法.前者是一种基于端元估计的解混方法, 代表性的研究有: Jia等[6]提出的添加分段光滑约束和稀疏约束的非负矩阵分解方法; Liu等[7]提出的基于$L_1-L_2$稀疏性约束的非负矩阵分解方法.此类方法的一个显著不足是需要进行端元估计; 而使用已知的光谱库作为端元集合的稀疏回归解混方法是目前研究的另一主要分支, 该方法的优点是不需要进行端元提取, 直接利用光谱库中光谱构成端元矩阵, 进而反演出丰度系数.由于光谱库中的光谱曲线个数远远大于实际场景中端元个数, 这样每一个观测到的混合像元光谱曲线在光谱库下的丰度系数向量就可以看成是"稀疏的".代表性研究有: 2011年Iordache[8]提出将光谱库引入线性高光谱解混模型中, 用已知的光谱库代替端元集合用于高光谱解混. Fu等[9]给出了OMP (Orthogonal matching pursuit)、ISMA (Iterative spectral mixture analysis)以及TwIST (Two step iterative shrinkage/thresholding algorithm)等算法进行稀疏性高光谱解的实验对比和应用分析, 但是这类方法的求解容易陷入局部最优, 从而影响解混的精度; 进而, Iordache等[10-11]提出一种基于交替迭代的SUnSAL (Sparse unmixing by variable splitting and augmented lagrangian)算法, 取得了较好的解混效果; 然而SUnSAL仅仅考虑了$l_1$稀疏性, 没有考虑空间信息, 为此2012年Iordache等[12-13]提出基于各向异性TV (Total variation)约束空间相关性的高光谱解混算法, 即SUnSAL-TV (Sparse unmixing by variable splitting and augmented lagrangian-total variation), 该方法能有效促进像元间的平滑性, 从而提高解混效果. 2013年宋义刚等[14]对基于光谱库的稀疏性高光谱解混方法和基于非负矩阵分解的稀疏性高光谱解混方法进行了综合分析和性能测试, 再次验证了稀疏性高光谱解混方法的有效性. 2014年Iordache等[15]提出了基于协同稀疏回归的高光谱解混算法CLSUnSAL (Collaborative SUnSAL), 利用丰度系数表现出的行稀疏特性, 施加协同稀疏, 提高解混效果.随后, Iordache等[16]又提出基于多重信号分类和协同稀疏回归的高光谱混合像元分解, 通过多信号分类算法产生比原始光谱库基数更小的端元子集, 进而利用协同稀疏性, 有效提高了解混效果. 2015年邓承志等[17]对于基于$L_1$稀疏正则化的高光谱混合像元分解算法比较, 表明全约束模型比非负约束模型能获得更精确的解混结果. 2016年Yang等[18]提出了一种耦合稀疏去噪和解混的CHyDU (Coupled HSI denoising and unmixing method)方法, 该方法利用低秩表示对丰度系数进行约束, 并且通过解混和去噪过程的相互促进来提高解混性能.虽然CHyDU方法的解混精度明显提高, 但是它却忽略了局部空间平滑性约束, 且模型求解相对复杂.与以上解混方法不同的另一种方法是Themelis等[19]在2012年提出的稀疏贝叶斯学习(Sparse Bayesian learning, SBL), 它基于贝叶斯思想, 对未知量进行概率建模, 并将丰度系数非负性加入先验模型中, 通过贝叶斯推断得到稀疏解, 但是该方法求解运算量大, 而且未考虑相邻像元中端元组合的联合稀疏性, 导致算法求解效率较低. 2016年孔繁锵等[20]提出一种基于复合正则化联合稀疏贝叶斯学习的高光谱稀疏解混算法CRMSBL (Compound regularized multiple sparse Bayesian learning), 该方法针对传统稀疏贝叶斯算法的不足, 考虑到稀疏丰度矩阵中非零元素所在位置具有结构化特征, 从而来提高解混精度.虽然CRMSBL方法的解混精度也明显提高, 但是它同样忽略了局部空间平滑性约束.

综上, 对于基于光谱库的稀疏回归解混方法, 有些算法仅考虑了高光谱图像(Hyperspectral image, HSI)的局部空间相关性, 如SUnSAL-TV方法, 它虽然考虑了高光谱图像的空间信息, 但是缺乏对丰度系数更为有效的稀疏性刻画; 有些算法仅考虑了协同稀疏性如CLSUnSAL方法, 该算法对丰度系数矩阵行稀疏约束, 有效刻画了丰度系数的行稀疏特性, 但是它并未考虑高光谱图像的局部空间平滑性, 导致局部空间信息的缺失, 影响解混精度.这两类算法的解混性能仍有待进一步提高, 从而本文将高光谱图像的局部空间相关性和全局协同稀疏性结合起来.既考虑高光谱图像相邻数据的局部相关性, 又利用丰度系数的行稀疏特性, 对丰度系数全局施加协同稀疏, 从而提高解混精度.由此, 本文提出一种基于协同稀疏和TV正则项的高光谱线性解混方法, 称为CLSUnSAL-TV (Collaborative SUnSAL-TV)方法, 该方法的流程图如图 1所示:首先, 建立基于光谱库的高光谱稀疏线性回归解混模型; 然后, 在丰度系数中加入协同稀疏和TV正则项, 用来刻画丰度系数的全局行稀疏性以及相邻像元的局部空间聚类特性, 提高稀疏解混的精确性; 最后, 采用交替方向乘子法[21] (Alternating direction method of multipliers, ADMM)求解模型, 得到解混结果.

图 1 基于协同稀疏和TV正则项的高光谱线性解混方法的流程图 Figure 1 Flowchart of the method based on collaborative sparsity and total variation
1 线性稀疏解混 1.1 线性光谱分解

线性混合模型假设在任何给定的光谱波段中, 每个像元光谱可由端元光谱线性组合而成[22], 即对于每个像元, 线性模型可以数学表示如下:

$ y_j=\sum\limits_{i=1}^q m_{ij}s_i+n_j $ (1)

其中, $y_j$是该混合像元在光谱波段$j$上的反射率, $m_{ij}$是端元$i$在光谱波段$j$的反射率, $s_i$是端元$i$对应混合像元的丰度系数值, $n_j$代表测量过程中产生的噪声, $q$代表端元数目.

假定高光谱遥感使用了$L$个光谱波段, 可以将式(1)写成如下的矩阵形式:

$ Y=MS+N $ (2)

其中, $Y\in{\bf{R}}^{L\times J}$表示高光谱图像($L$为波段数, $J$为每个波段图像中的像素个数); $M=[\pmb{m}_1, \pmb{m}_2, \cdots, \pmb{m}_q]\in{\bf R}^{L\times q}$表示端元矩阵, 代表第$i$个端元; $S=[\pmb{s}_1, \pmb{s}_2, \cdots, \pmb{s}_J]^\mathit{\boldsymbol{T}}\in{\bf R}^{q\times J}$表示丰度系数矩阵, $\pmb{s}_j\in{\bf R}^{q\times 1}$代表像元$j$中各个端元所占比例系数; $N\in{\bf R}^{L\times L}$是系统噪声.解混的目的是从观测的高光谱数据$Y\equiv{\{\pmb{y}_j\in{\bf R}^{L\times 1}, j=1, \cdots, J\}}$中估计出端元矩阵$M$及其对应的丰度系数矩阵$S$.

根据丰度的物理意义, 丰度系数[23]满足"非负性"约束(Abundance non-negativity constraint, ANC): $S_{ij}\geq0, \forall i=1, 2, \cdots , q, \forall j=1, 2, \cdots, J$, 并且满足"和为一"约束(Abundance sum-to-one constraint, ASC): $\mathit{1}^\mathit{\boldsymbol{T}}\pmb{s}_j =1, \forall j=1, 2, \cdots, J$, 其中$\mathit{1}^\mathit{\boldsymbol{T}}\in{\bf R}^{L\times 1}$且各元素全为1.

1.2 稀疏性分解

假设在式(2)中所提出的稀疏分解中使用一个先验的光谱特征库, 即用已知的光谱库代替端元矩阵, 则线性混合模型可以改写为

$ Y=AX+N $ (3)

其中, $A\in{\bf{R}}^{L\times g}$, $g$为光谱库中地物光谱数目, $X\in{\bf{R}}^{g\times J}$为光谱库$A$对应的丰度系数矩阵.

根据基于光谱库的高光谱混合像元分解机理, $A$中的谱线数远大于式(2)中$M$的端元个数, 即$X$中非零系数的个数$q\leq g$, 通常$X$的每一列是一个$k$稀疏向量, 即它有$k$个非零元素.这样形成的稀疏性能够提高高光谱混合像元解混的准确性, 由此得到约束的$l_0$优化问题:

$ \begin{array}{l l } \min\limits_{X}&\dfrac{1}{2}\|AX-Y\|_{\text{F}}^{2}+\lambda\|X\|_{0}\\ &\text{s}.\ \text{t}.~X\geq0, \mathit{1}^{\mathit{\boldsymbol{T}}}X=\mathit{1} \end{array} $ (4)

其中, 第一项为重构误差项, 第二项为稀疏性约束项; $\|X\|_{\text{F}}=\sqrt{{\rm tr}{\{XX^{\mathit{\boldsymbol{T}}}\}}}$为Frobenius范数; $\|X\|_{0}$代表了$X$$l_0$范数, 它计算了$X$中非零元素的个数; $X\geq0$, 即矩阵$X$中每个元素均非负; $\mathit{1}^{\mathit{\boldsymbol{T}}}X=\mathit{1}$代表$X$的每一列满足"和为一''约束, $\mathit{1}^\mathit{\boldsymbol{T}}\in{\bf R}^{L\times 1}$且各元素全为1; $\lambda$是一个平衡保真项和正则项的参数.

根据稀疏表示理论, $l_0$范数的最小化是NP难且非凸问题, 求解困难.因此常用$l_1$范数替代$l_0$范数, 且当矩阵$A$满足RIP (Restricted isometry property)条件时[8, 10], $l_0$问题等价于$l_1$凸优化问题:

$ \begin{array}{l l } \min\limits_{X}&\dfrac{1}{2}\|AX-Y\|_{\text{F}}^{2}+\lambda\|X\|_{1, 1}\\ &\text{s}.\ \text{t}.~X\geq0, \mathit{1}^{\mathit{\boldsymbol{T}}}X=\mathit{1} \end{array} $ (5)

其中, $\|X\|_{1, 1}=\sum_{k=1}^g\|\pmb{x}^k\|_{1}$ ($\pmb{x}^k$代表$X$的第$k$行)表示矩阵元素之和.

2 本文提出的稀疏解混方法 2.1 联合协同稀疏和TV正则项的线性解混模型

高光谱图像局部相邻像元通常包含相似的物质, 其光谱特征具有相似性, 利用高光谱图像像元间的这种相关性, 可以有效提高解混的准确性.进一步分析像元间的空间相关性, 可以发现这种相关性通常表现为像元间的平滑性, 由此加入TV正则项的稀疏回归模型[12]可以表示如下:

$ \begin{array}{l l l} \min\limits_{X}&\dfrac{1}{2}\|AX-Y\|_{\text{F}}^{2}+\lambda\|X\|_{1, 1}+\lambda_{\text{TV}}\text{TV}(X)\\ &{\rm s.\, t.}~X\geq0, \mathit{1}^{\mathit{\boldsymbol{T}}}X=\mathit{1} \end{array} $ (6)

其中, $\text{TV}(X)\equiv\sum_{(i, j)\in \varepsilon}\|\pmb{x}_{i}-\pmb{x}_{j}\|_{1}$是各向异性TV, 它刻画了相邻像元间的平滑性约束, 且$\varepsilon$表示图像中水平和垂直相邻像元的集合, $\lambda$$\lambda_{\text{TV}}$为正则化参数.

以上模型的稀疏性是添加在丰度系数上的, 由于光谱库足够大, 而场景中的端元个数远远小于给定光谱库中光谱个数, 所以丰度系数必然从全局上表现出行稀疏的特性.行稀疏性比$l_1$稀疏更能有效刻画丰度系数的稀疏模式, 具有更强的稀疏性, 同时, 耦合局部空间平滑性约束, 有效提升解混能力.基于协同稀疏和TV正则化相结合的解混模型可表示如下:

$ \begin{align} \min\limits_{X}&\dfrac{1}{2}\|AX-Y\|_{\text{F}}^{2}+ \lambda\|X\|_{2, 1}+\lambda_{\text{TV}}\text{TV}(X)\nonumber\\ &{\rm s.\, t.}~X\geq0, \mathit{1}^{\mathit{\boldsymbol{T}}}X=\mathit{1} \end{align} $ (7)

其中, $\|X\|_{2, 1}=\sum_{k=1}^g\|\pmb{x}^k\|_{2}$$l_{2, 1}$混合范数, 它刻画了$X$的行稀疏性.

定义$H_{h}$: ${\bf R}^{g\times J}\rightarrow {\bf R}^{g\times J }$$X$中像元和其相邻像元间水平差分的线性算子[12], 即$H_{h}X=[{\pmb d}_{1}, {\pmb d}_{2}, \cdots, {\pmb d}_{n}]$, ${\pmb d}_{i}= {\pmb x}_{i}-{\pmb x}_{ih}$, $i$$ih$分别代表像元$i$和它的水平相邻像元.同理, 定义一个垂直线性差分算子$H_{v}$: ${\bf R}^{g\times J}\rightarrow {\bf R}^{g\times J }$, 即$H_{v}X=[{\pmb v}_{1}, {\pmb v}_{2}, \cdots, {\pmb v}_{n}]$, ${\pmb v}_{i}= {\pmb x}_{i}-{\pmb x}_{iv}$, $i$$iv$分别代表像元$i$和它的垂直相邻像元.基于以上定义, 我们有$HX\equiv\left( \begin{array}{c} H_{h}X\\ H_{v}X\\ \end{array} \right)$.进而, 式(7)可重写为

$ \begin{align} \min\limits_{X}&\dfrac{1}{2}\|AX-Y\|_{\text{F}}^{2}+\lambda\|X\|_{2, 1}+\nonumber\\ &\lambda_{\text{TV}}\parallel HX\parallel_{1, 1}+l_{\bf{R}+}(X)+l_{1}(\mathit{1}^{\mathit{\boldsymbol{T}}}X) \end{align} $ (8)

其中, $l_{\bf{R}+}(X)$表示"非负性"约束的示性函数(其中$l_{\bf{R}+}(\cdot)$定义为:当$a\geq 0$时, $l_{\bf{R}+}(a)=0$, 当$a<0$$l_{\bf{R}+}(a)=\infty$). $l_{1}(\mathit{1}^{\mathit{\boldsymbol{T}}}X)$表示"和为一''约束的示性函数. $l_{\bf{R}+}(X)$, $l_{\bf{R}+}(a)=\infty$). $l_{1}(\mathit{1}^{\mathit{\boldsymbol{T}}}X)$分别保证了解混过程中丰度系数$X$的"非负性"约束和"和为一"约束得以满足. $\lambda\geq0$$\lambda_{\text{TV}}\geq0$分别为$l_{2, 1}$稀疏约束项和TV正则项的正则化参数, 两者平衡了$l_{2, 1}$稀疏约束项、TV正则项与其他各项之间的关系.当$\lambda_{\text{TV}}=0$时, 模型(8)即为经典的CLSUnSAL方法.理论上随着噪声的增强, $\lambda$$\lambda_{\text{TV}}$取值相对较大时才能获得最优的解, 且随着的$\lambda\geq0$$\lambda_{\text{TV}}\geq0$增大, 其对应正则项的作用在解混中增强.

2.2 模型求解

对于式(8)中优化问题, 本文采用交替方向乘子法(ADMM)求解, 通过引入变量$V_1, V_2, V_3, V_4, V_5, V_6$, 可以得到如下等价形式:

$ \begin{array}{l l } \min\limits_{X, V_1, V_2, V_3, V_4, V_5, V_6}&\dfrac{1}{2}\parallel V_1-Y\parallel_{\text{F}}^2+{\lambda}\parallel V_2\parallel_{2, 1}+\nonumber\\ &\lambda_{\text{TV}}\parallel V_{4}\parallel_{1}+l_{\bf{R}+}(V_5)+ \\ &l_{1}(\mathit{1}^{\mathit{\boldsymbol{T}}}V_6)\end{array} $
$ \begin{array}{l l } &V_1=AX\\ &V_2=X\\ &V_3=X\\ \text{s}.\ \text{t}.&V_4=HV_3\\ &V_5=X\\ &V_6=X \end{array} $ (9)

式(9)中将$V_3=HX$替换为$V_3=X$$V_4=HV_3$, 这样可以先求解出$V_3$, 再通过$V_3$求解出$V_4$, 降低原问题的复杂度.

最优化问题(9)是一个典型的约束最小化问题, 其增广拉格朗日方程为

$ \begin{align} \ell(X, &V_1, V_2, V_3, V_4, V_5, V_6, D_1, D_2, D_3, D_4, D_5, D_6)\equiv\nonumber\\ &\dfrac{1}{2}\parallel V_1-Y\parallel_{\text{F}}^2+{\lambda}\parallel V_2\parallel_{2, 1}+\lambda_{\text{TV}}\parallel V_{4}\parallel_{1, 1}+\nonumber\\ &l_{\bf{R}+}(V_5)+l_1(\mathit{1}^{\mathit{\boldsymbol{T}}}V_6)+\dfrac{\mu}{2}\parallel AX-V_1-D_1\parallel_\text{F}^{2}+\nonumber\\ &\dfrac{\mu}{2}\parallel X-V_2-D_2\parallel _\text{F}^{2}+ \dfrac{\mu}{2}\parallel X-V_3-D_3\parallel_\text{F}^{2}+\nonumber\\ &\dfrac{\mu}{2}\parallel HV_3-V_4-D_4\parallel_\text{F}^{2} +\nonumber\\&\dfrac{\mu}{2}\parallel X-V_5-D_5\parallel _\text{F}^{2}+ \dfrac{\mu}{2}\parallel X-V_6-D_6\parallel _\text{F}^{2} \end{align} $ (10)

其中, $\mu$为一个常数, $D_1, D_2, D_3, D_4, D_5, D_6$分别表示了6个拉格朗日乘子.具体算法如下:

步骤1. 设置$k=0$, 选择$\mu>0$, 初始化$X^{(0)}$, $V^{(0)}_{1}, \cdots, V^{(0)}_{6}$, $D^{(0)}_{1}, \cdots, D^{(0)}_{6}$.

步骤2. 计算$X^{(k+1)}\leftarrow\arg\min_{X}\ell(X, V_{1}^{(k)}, \cdots, \\V_{6}^{(k)}, D_{1}^{(k)}, \cdots, D_{6}^{(k)})$.

步骤3. 计算$V_{i}^{(k+1)}\leftarrow\arg\min_{V_{i}}\ell(X^{(k+1)}, V_{1}^{(k)}, \\\cdots, V_{i}^{(k)}, \cdots, V_{6}^{(k)}), i=1, \cdots, 6$.

步骤4. 更新拉格朗日乘子

$D_{1}^{(k+1)}\leftarrow D_{1}^{(k)}-AX^{(k+1)}+V_{1}^{(k+1)}$;

$D_{4}^{k+1}\leftarrow D_{4}^{(k)}-HV_{3}^{(k+1)}+V_{4}^{(k+1)}$;

$D_{i}^{(k+1)}\leftarrow D_{i}^{(k)}-X^{(k+1)}+V_{i}^{(k+1)}, i=2, 3, 5, 6$.

步骤5. 如果满足停止条件, 结束算法; 否则, $k\leftarrow k+1$, 算法转至步骤2.

对于上述算法, 步骤2中每一步的$X$值可以由如下式子得到:

$ \begin{align} X^{(k+1)}&\leftarrow\arg \min\limits_{X}\dfrac{\mu}{2}\parallel AX-V_1^{(k)}-D_1^{(k)}\parallel^{2}_{\text{F}}+\nonumber\\ &\dfrac{\mu}{2}\parallel X-V_i^{(k)}-D_i^{(k)}\parallel _\text{F}^{2}, ~~i=2, 3, 5, 6\end{align} $ (11)

通过对式(11)求导可以得到$X$的解:

$ \begin{align} X^{(k+1)}&\leftarrow(A^{\mathit{\boldsymbol{T}}}A+4I)^{-1} (A^{\mathit{\boldsymbol{T}}}\boldsymbol\xi_{1}+\boldsymbol\xi_{2}+ \boldsymbol\xi_{3}+\nonumber\\ &\boldsymbol\xi_{5}+\boldsymbol\xi_{6}) \end{align} $ (12)

其中, $I$为单位矩阵, $A^{\mathit{\boldsymbol{T}}}$$A$的转置, $\boldsymbol\xi_{i}=V_{i}^{(k)}+D_{i}^{(k)}, i=1, 2, 3, 5, 6$.

该算法的步骤3中计算了每个变量$V_1, V_2, V_3, V_4, V_5, V_6$的值.为了计算$V_1$, 优化问题可以写为

$ \begin{align} V_1^{(k+1)}&\leftarrow\arg \min\limits_{V_{1}} \dfrac{1}{2}\parallel V_1-Y\parallel _\text{F}^{2}+\nonumber\\&\dfrac{\mu}{2}\parallel AX^{(k+1)}-V_1-D_1^{(k)}\parallel^{2}_{\text{F}} \end{align} $ (13)

其解为: $V_1^{(k+1)}\leftarrow\frac{1}{1+\mu}[Y+\mu(AX^{(k+1)}-D_{1}^{(k)})]$.

为了计算$V_2$, 优化问题可以写为

$ \begin{align} V_2^{(k+1)}&\leftarrow\arg \min\limits_{V_{2}} {\lambda}\parallel V_2\parallel _{2, 1}+\nonumber\\ &\dfrac{\mu}{2}\parallel X^{(k+1)}-V_2-D_2^{(k)}\parallel^{2}_{\text{F}} \end{align} $ (14)

其解可由行软阈值Vect-soft法[15]求得: $V_2^{(k+1)}\leftarrow \mbox{vect-soft}(\varphi_{2}, {\lambda}/{\mu})$, 其中$\varphi_{2}=X^{(k+1)}-D_{2}^{(k)}$$\mbox{vect-soft}(\varphi_{2}, {\lambda}/{\mu})=\varphi_{2}(\frac{\max\{\|\varphi_{2}\|_{2}-\frac{\lambda}{\mu}, 0\}}{\max\{\|\varphi_{2}\|_{2}-\frac{\lambda}{\mu}, 0\}+\frac{\lambda}{\mu}})$.

为了计算$V_3$, 优化问题可以写为

$ \begin{align} V_3^{(k+1)}&\leftarrow\arg \min\limits_{V_{3}} \dfrac{\mu}{2}\parallel X^{(k+1)}-V_3-D_{3}^{(k)}\parallel _\text{F}^{2}+\nonumber\\&\dfrac{\mu}{2}\parallel HV_{3}-V_4^{(k)}-D_4^{(k)}\parallel^{2}_{\text{F}} \end{align} $ (15)

其解为: $V_3^{(k+1)}\leftarrow(H^{\mathit{\boldsymbol{T}}}H+I)^{-1}(X^{(k+1)}-D_{3}^{(k)}+H^{\mathit{\boldsymbol{T}}}\boldsymbol\xi_{4})$, 其中$\boldsymbol\xi_{4}=V_{4}^{(k)}+D_{4}^{(k)}$, 通过傅里叶快速算法求解得出$V_{3}^{(k+1)}\leftarrow{\cal F}^{-1}(\frac{{\cal F}(X^{(k+1)}-D_{3}^{(k)}+H^{\mathit{\boldsymbol{T}}}\boldsymbol\xi_{4})}{(1+{\cal F}(H_{h})^{\mathit{\boldsymbol{T}}}{\cal F}(H_{h})+{\cal F}(H_{v})^{\mathit{\boldsymbol{T}}}{\cal F}(H_{v})})$, 其中${\cal F}(\cdot)$为傅里叶变换, ${\cal F}^{-1}(\cdot)$为傅里叶逆变换.

$V_{4}$的计算可以通过解决如下优化问题得到:

$ \begin{align} V_4^{(k+1)}&\leftarrow\arg \min\limits_{V_{4}} \lambda_{\text{TV}}\parallel V_4\parallel _{1, 1}+\nonumber\\&\dfrac{\mu}{2}\parallel HV_{3}^{(k)}-V_4-D_4^{(k)}\parallel^{2}_{\text{F}} \end{align} $ (16)

其解可由软阈值Soft-threshold法[20]求得: $V_4^{(k+1)}\leftarrow \mathit{\boldsymbol{soft}}(D_{4}^{(k)}-HV_{3}^{(k)}, {\lambda_{\text{TV}}}/{\mu})$, 其中$\mathit{\boldsymbol{soft}}(D_{4}^{(k)}-HV_{3}^{(k)}, {\lambda_{\text{TV}}}/{\mu})$$\mathit{\boldsymbol{sgn}}(D_{4}^{(k)}-HV_{3}^{(k)})\cdot\max(|D_{4}^{(k)}-HV_{3}^{(k)}|-{\lambda_{\text{TV}}}/{\mu}, 0)$.

对于$V_{5}$的优化问题:

$ \begin{align} V_5^{(k+1)}&\leftarrow\arg \min\limits_{V_{5}} l_{\bf{R}^{+}}(V_5)+\nonumber\\ &\dfrac{\mu}{2}\| X^{(k+1)}-V_5-D_5^{(k)}\|^{2}_{\text{F}} \end{align} $ (17)

其中, $V_{5}$的值可以由$V_{5}^{(k+1)}\leftarrow \max(X^{(k+1)}-D_{5}^{(k)}, 0)$得到.

最后, 对于$V_{6}$的优化问题:

$ \begin{align} V_6^{(k+1)}&\leftarrow\arg \min\limits_{V_{6}} l_{1} (\mathit{1}^{\mathit{\boldsymbol{T}}}V_{6})+\nonumber\\ &\dfrac{\mu}{2}\| X^{(k+1)}-V_6-D_6^{(k)}\|^{2}_{\text{F}} \end{align} $ (18)

其解为$V_6^{(k+1)}\leftarrow(X^{(k+1)}-D_6^{(k)})+(Q\otimes W)$, 其中$W=\frac{1}{n}[1-\sum_{i=1}^{n}(X^{(k+1)}-D_6^{(k)})]$, $Q=\mathit{\boldsymbol{ones}}(n, 1)$, $\otimes$为克罗内克积(Kronecker product).

在上述整个算法过程中, 最主要的步骤是$X$$V_3$的求解, 两者的复杂度分别为O$ (L^{2}J)$和O$ (LJ\log J)$, 其中$L$为光谱波段数目, $J$为图像中像元数目.因此, 整个算法每次迭代的复杂度为O$ (LJ\cdot\max\{L, \log J\})$.

3 实验与分析 3.1 模拟数据实验

在本节中, 我们采用模拟数据测试算法性能, 本文提出的算法(CLSUnSAL-TV)将与基于$l_1$范数的两种算法[10] (式(5)中$\lambda=0$时的算法NCLS, $\lambda\neq0$时的算法SUnSAL), 基于协同稀疏的算法[15] (CLSUnSAL), 以及基于TV正则化的两种算法[12] (式(6)中$\lambda=0$时的算法NCLS-TV, $\lambda\neq0$时的算法SUnSAL-TV)进行比较.解混精度采用SRE信号重建误差[8, 10, 13]以及均方根误差RMSE[24]来度量. SRE数值越大, RMSE数值越小, 说明算法的解混精度越高, 其定义分别如下式(19)和式(20):

$ \text{SRE(dB)}\equiv 10\lg\frac{{\text{E}}[\|X\|_{2}^{2}]}{{\text{E}}[\|X-\hat{X}\|_{2}^{2}]} $ (19)
$ \text{RMSE}\equiv\sqrt{\frac{1}{g\times J}\sum\limits_{i=1}^{g}\sum\limits_{j=1}^{J}(X_{ij}-\widehat{X}_{ij})^{2}} $ (20)

其中, $X$为真实丰度系数, $\hat{X}$为解混算法重建的丰度系数; $g$为光谱库中光谱数, $J$为像元总数, $X_{ij}$为真实丰度系数矩阵中的元素, $\widehat{X}_{ij}$为解混算法重建的丰度系数矩阵中的元素.

本节的测试环境为Intel Core i5 CPU 3.20 GHz、内存4.0 GB、Matlab R2013b.本节算法的参数依据经验设置为如下: $\lambda, \lambda_{\text{TV}}\in$ [0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5], 且$\mu\in$ [0.01, 0.02, 0.03, 0.04, 0.05, 0.08, 0.1, 0.2, 0.3, 0.4, 0.5, 0.8, 1].从这些范围中选择合适的$\lambda$, $\lambda_{\text{TV}}$以及$\mu$, 得到最佳的SRE值和RMSE值, 并用粗体显示.

3.1.1 模拟数据集

实验中采用的光谱库为$A\in {\bf{R}}^{224\times240}$: $A$由在2007年9月发布的USGS光谱库中选择的240种地物光谱构成, 光谱范围为0.4 $\sim 2.5 \rm\mu$m.

从光谱库$A$中随机选取5个光谱特征作为端元线性生成模拟数据集, 数据集有$75\, \times\, 75$个像元且每个像元有224个波段, 图 2中显示了5种端元的光谱曲线图. 图 3中分别显示了模拟数据图像以及5个端元的真实丰度图, 模拟数据图像中既有纯像元, 又有由2至5个端元混合成的混合像元, 这些像元均匀分布在固定的方形区域中.模拟数据的背景像元由5个相同的端元混合组成, 但其各自的丰度系数固定为0.1149, 0.0742, 0.2003, 0.2055, 0.4051.按照上述的描述过程产生数据集后, 在高斯白噪声污染的情况下进行实验.实验中采用3种不同信噪比的噪声, (信噪比$\text{SNR}\equiv\frac{{\text{E}}\|AX\|^{2}_{2}}{{\text{E}}\|N\|^{2}_{2}}$)即20 dB、30 dB、40 dB.

图 2 5个端元在光谱库中的谱线 Figure 2 Spectral characteristic curves of five endmembers
图 3 模拟数据图以及其中5个端元的真实丰度图 Figure 3 Simulated image and five true fractional abundances of endmembers in the simulated dataset
3.1.2 实验结果与分析

在白噪声情况下, 我们采取6种稀疏解混方法在光谱库$A$中对模拟数据集进行解混, 通过计算其SRE和RMSE来比较6种方法的解混精度. 表 1中显示了在添加"和为一"约束和不添加"和为一"约束时, 6种方法解混得到的SRE和RMSE结果以及解混时间$t\, (s), $同时图 4中给出了不同噪声下CLSUnSAL-TV算法得到的关于参数$\lambda$$\lambda_{\text{TV}}$的SRE函数图.由表 1结合图 4可以看出: $\lambda$$\lambda_{\text{TV}}$的值对解混的结果具有较大影响, 通常噪声较大时, 最优结果对应的$\lambda$$\lambda_{\text{TV}}$的取值相对较大:相同噪声情况下, 添加"和为一"约束的SRE值均高于不添加"和为一"约束的SRE值, 且添加"和为一"约束的RMSE值均不高于不添加"和为一"约束的RMSE值, 从而表明添加"和为一"约束的解混效果比不添加"和为一"约束的解混效果好; NCLS-TVSUnSAL-TV、CLSUnSAL-TV的SRE值远高于其他3种算法, 且NCLS-TV、SUnSAL-TV、CLSUnSAL-TV的RMSE值远低于其他3种算法, 这意味着TV正则项有效刻画了地物的局部空间聚类性; 在不同噪声情况下, CLSunSAL-TV算法的SRE值最大, RMSE值最小, 表明全局行稀疏性进一步提高解混能力, 使得解混效果最佳; 在不同噪声下, 与当前较好的解混方法SUnSAL-TV相比, CLSunSAL-TV算法的运行时间均小于SUnSAL-TV算法的运行时间, 这表明CLSunSAL-TV算法在时间上有一定的优越性.

表 1 6种方法在最优参数($\lambda, \lambda_{\text{TV}}, \mu$)时对模拟数据集解混得到的SRE (dB), RMSE ($10^{-2}$)以及所需时间$t$ (s) Table 1 SRE, RMSE ($10^{-2}$), the running time (seconds), and optimal parameters of the six unmixing methods for the simulated dataset
图 4 不同噪声下CLSUnSAL-TV算法得到的关于参数$\lambda$$\lambda_{\text{TV}}$的SRE函数图 Figure 4 SRE (dB) as a function of parameters $\lambda$ and $\lambda_{\text{TV}}$ obtained by CLSUnSAL-TV algorithm in different noise levels

为了更直观地显示解混的效果, 图 5显示了$\text{SNR}=40 \mathit{\boldsymbol{dB}}$时各种方法解混得到的丰度图.这里我们仅仅给出3种端元(端元1、端元2、端元5)的估计丰度图, 因为另外两种端元表现了相似的解混效果.由于添加"和为一"约束的解混效果比不添加"和为一"约束的解混效果好, 所以在图 5显示提取的端元丰度图时, 仅考虑添加"和为一"约束的情况. 图 5中的第一行显示了3种端元在光谱库中的谱线图, 这3种端元分别为Jarosite GDS99、Jarosite GDS101、Alunite 3种矿物, 图 5中的其余行分别给出6种算法在信噪比SNR为40 dB时重建的丰度图.从图 5看出:在NCLS、SUnSAL、CLSUnSAL解混的丰度图中有很多噪声点, 而施加了TV正则项的NCLS-TV、SUnSAL-TV、CLSUnSAL-TV 3种方法得到的解混图中几乎没有噪声点, 这充分表明TV正则项对噪声抑制较好, 有利于提高解混精度; CLSUnSAL-TV算法得到的丰度图中不但几乎没有噪声, 而且比其他方法有着更多小正方形区域的丰度信息, 丰度图像最接近真实丰度图, 因此CLSUnSAL-TV方法的混合像元分解精度最高.

图 5 信噪比为40 dB情况下6种算法解混得到的丰度图 Figure 5 Reconstructed abundances of six methods for simulated dataset when SNR = 40 (dB)
3.2 真实数据实验 3.2.1 Cuprite数据集

实际数据采用AVIRIS光谱仪采集的Cuprite数据, 该数据集为1995年7月获取的美国内华达州Cuprite采矿区数据.图像大小为$250\, \times\, 191$, 该场景包含0.4到2.5 $\mu$m之间的224个光谱波段, 其光谱分辨率为10 nm.去除水汽干扰和噪声影响的波段1 $\sim$ 2, 105 $\sim$ 115, 150 $\sim$ 170, 223 $\sim$ 224, 剩余188个有效波段.在图 6中显示了Cuprite数据的RGB伪彩合成图.

图 6 Cuprite数据的RGB伪彩合成图(波段15、23、117) Figure 6 The false color image of AVIRIS Cuprite (band 15、23、117)
3.2.2 实验结果与分析

基于USGS光谱库$A$, 对Cuprite数据使用6种稀疏解混方法得到端元的丰度系数图. 图 7中显示了利用以上6种稀疏分解算法提取3种最具代表性矿物的丰度图, 这3种矿物分别为Alunite (明矾石)、Buddingtonite (水铵长石)、Chalcedony (玉髓). 图 7中第1行显示了这3种矿物在光谱库中的谱线图, 第2行给出Cuprite数据集中3种矿物的实际丰度图, 用作参照对比. 图 7中第3行到第8行分别显示了应用NCLS、SUnSAL、CLSUnSAL、NCLS-TV、SUnSAL-TV、CLSUnSAL-TV 6种方法解混得到的矿物丰度估计图.根据图 7总结出:从第3行到第8行中的矿物丰度估计图越来越接近真实的矿物丰度图, 即解混效果越来越好: NCLS-TV、SUnSAL-TV、CLSUnSAL-TV表现了相似的矿物丰度估计图, 它们与NCLS、SUnSAL、CLSUnSAL方法相比, 丰度估计图中有着更少的离群值, 丰度系数间更加平滑; 尽管很难对SUnSAL-TV和CLSUnSAL-TV算法解混的矿物丰度图做出比较, 但CLSUnSAL-TV解混得到的Alunite (明矾石)和Buddingtonite (水铵长石)的丰度图比SUnSAL-TV解混得到的更接近于参照图, 即表明CLSUnSAL-TV方法的解混精度最高, 体现了本文算法的优越性.

图 7 对真实数据使用6种方法解混得到的丰度图 Figure 7 Reconstructed abundances of six methods for real dataset
4 结论

同时考虑丰度系数的全局行稀疏性和相邻像元之间的局部空间平滑特性, 本文提出了一种基于协同稀疏和TV正则项的高光谱空谱联合解混方法, 即CLSUnSAL-TV.该方法与CLSUnSAL和SUnSAL-TV相比, 具有更佳的解混效果, 表明同时考虑丰度系数的全局行稀疏性和相邻像元之间的局部空间平滑性对基于已知光谱库的稀疏回归解混方法的必要性, 克服了以往对于解混模型只重视挖掘空间相关性而忽视稀疏性精确刻画的缺点, 为以后的空谱联合稀疏性回归解混方法提供了新思路.模拟和真实高光谱数据实验结果验证了本文方法是精确且有效的.

参考文献
1
Tong Qing-Xi, Zhang Bing, Zheng Lan-Feng. Hyperspectral Remote Sensing:Theory, Technology and Applications. Beijing: Higher Education Press, 2006.
( 童庆禧, 张兵, 郑兰芬. 高光谱遥感:原理、技术与应用. 北京: 高等教育出版社, 2006.)
2
Yang Hua-Dong. Research on spectral unmixing algorithms for hyperspectral remote sensing image[Ph. D. dissertation], Dalian Maritime University, China, 2015
(杨华东. 高光谱遥感影像光谱解混算法研究[博士学位论文], 大连海事大学, 中国, 2015) http://cdmd.cnki.com.cn/Article/CDMD-10151-1015354106.htm
3
Bioucas-Dias J M, Nascimento J M P. Hyperspectral subspace identification. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(8): 2435-2445. DOI:10.1109/TGRS.2008.918089
4
Bioucas-Dias J M, Plaza A, Dobigeon N, Parente M, Du Q, Gader P, Chanussot J. Hyperspectral unmixing overview:geometrical, statistical, and sparse regression-based approaches. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012, 5(2): 354-379. DOI:10.1109/JSTARS.2012.2194696
5
Wang Mao-Zhi, Xu Wen-Xi, Wang Lu, Guo Ke. Research progress on endmember extraction algorithm and its classification of hyperspectral remote sensing imagery. Remote Sensing Technology and Application, 2015, 30(4): 616-625.
( 王茂芝, 徐文皙, 王璐, 郭科. 高光谱遥感影像端元提取算法研究进展及分类. 遥感技术与应用, 2015, 30(4): 616-625.)
6
Jia S, Qian Y T. Constrained nonnegative matrix factorization for hyperspectral unmixing. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(1): 161-173. DOI:10.1109/TGRS.2008.2002882
7
Liu J J, Wu Z B, Wei Z H, Xiao L, Sun L. A novel sparsity constrained nonnegative matrix factorization for hyperspectral unmixing. In: Proceedings of the 2012 IEEE International Conference on Geoscience and Remote Sensing Symposium (IGARSS). Munich, Germany: IEEE, 2012. 1389-1392
8
Iordache M D. A sparse regression approach to hyperspectral unmixing[Ph. D. dissertation], Universidade Tecnica de Lisboa, Portugal, 2011
9
Fu X, Ma W K, Chan T H, Bioucas-Dias J M. Self-dictionary sparse regression for hyperspectral unmixing:greedy pursuit and pure pixel search are related. IEEE Journal of Selected Topics in Signal Processing, 2015, 9(6): 1128-1141. DOI:10.1109/JSTSP.2015.2410763
10
Iordache M D, Bioucas-Dias J M, Plaza A. Sparse unmixing of hyperspectral data. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(6): 2014-2039. DOI:10.1109/TGRS.2010.2098413
11
Bioucas-Dias J M. A variable splitting augmented lagrangian approach to linear spectral unmixing. In: Proceedings of the 1st IEEE Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing. Grenoble, France: IEEE, 2009. 1-4
12
Iordache M D, Bioucas-Dias J M, Plaza A. Total variation spatial regularization for sparse hyperspectral unmixing. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(11): 4484-4502. DOI:10.1109/TGRS.2012.2191590
13
Sun L, Wu Z B, Liu J J, Xiao L, Wei Z H. Supervised spectral-spatial hyperspectral image classification with weighted Markov random fields. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(3): 1490-1503. DOI:10.1109/TGRS.2014.2344442
14
Song Yi-Gang, Wu Ze-Bin, Wei Zhi-Hui, Sun Le, Liu Jian-Jun. Survey of sparsity constrained hyperspectral unmixing. Journal of Nanjing University of Science and Technology, 2013, 37(4): 486-492.
( 宋义刚, 吴泽彬, 韦志辉, 孙乐, 刘建军. 稀疏性高光谱解混方法研究. 南京理工大学学报, 2013, 37(4): 486-492.)
15
Iordache M D, Bioucas-Dias J M, Plaza A. Collaborative sparse regression for hyperspectral unmixing. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(1): 341-354. DOI:10.1109/TGRS.2013.2240001
16
Iordache M D, Bioucas-Dias J M, Plaza A, Somers B. MUSIC-CSR:hyperspectral unmixing via multiple signal classification and collaborative sparse regression. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(7): 4364-4382. DOI:10.1109/TGRS.2013.2281589
17
Deng Cheng-Zhi, Zhang Shao-Quan, Wang Sheng-Qian, Tian Wei, Zhu Hua-Sheng, Hu Sai-Feng. Hyperspectral unmixing algorithm based on L1 regularization. Infrared and Laser Engineering, 2015, 44(3): 1092-1097.
( 邓承志, 张绍泉, 汪胜前, 田伟, 朱华生, 胡赛凤. L1稀疏正则化的高光谱混合像元分解算法比较. 红外与激光工程, 2015, 44(3): 1092-1097.)
18
Yang J X, Zhao Y Q, Chan J C W, Kong S G. Coupled sparse denoising and unmixing with low-rank constraint for hyperspectral image. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(3): 1818-1833. DOI:10.1109/TGRS.2015.2489218
19
Themelis K E, Rontogiannis A, Koutroumbas K D. A novel hierarchical Bayesian approach for sparse semisupervised hyperspectral unmixing. IEEE Transactions on Signal Processing, 2012, 60(2): 585-599. DOI:10.1109/TSP.2011.2174052
20
Kong Fan-Qiang, Guo Wen-Jun, Shen Qiu, Wang Dan-Dan. Compound regularized multiple sparse Bayesian learning algorithm for sparse unmixing of hyperspectral data. Waves Journal of Infrared and Millimeter Waves, 2016, 35(2): 219-226.
( 孔繁锵, 郭文骏, 沈秋, 王丹丹. 复合正则化联合稀疏贝叶斯学习的高光谱稀疏解混算法. 红外与毫米波学报, 2016, 35(2): 219-226. DOI:10.11972/j.issn.1001-9014.2016.02.018)
21
Bioucas-Dias J M, Figueiredo M A T. Alternating direction algorithms for constrained sparse regression: application to hyperspectral unmixing. In: Proceedings of the 2nd Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS). Reykjavik, Iceland: IEEE, 2010. 1-4
22
Altmann Y, Pereyra M, Bioucas-Dias J M. Collaborative sparse regression using spatially correlated supports-application to hyperspectral unmixing. IEEE Transactions on Image Processing, 2015, 24(12): 5800-5811. DOI:10.1109/TIP.2015.2487862
23
Zhao Chun-Hui, Xiao Jian-Yu, Qi Bin. An improved OMP algorithm for hyperspectral sparse unmixing. Journal of Shenyang University (Natural Science), 2015, 27(3): 206-213.
( 赵春晖, 肖健钰, 齐滨. 一种改进的OMP高光谱稀疏解混算法. 沈阳大学学报(自然科学版), 2015, 27(3): 206-213.)
24
Guerra R, Santos L, López S, Sarmiento R. A new fast algorithm for linearly unmixing hyperspectral images. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(12): 6752-6765. DOI:10.1109/TGRS.2015.2447573