高光谱图像非负稀疏分量分解建模与鲁棒性解混方法

汪顺清，杨劲翔，邵远天，肖亮

南京理工大学计算机科学与工程学院,南京 210094

0 引 言

高光谱遥感数据包含丰富的空间维和光谱维信息，具有图谱合一、光谱分辨率高和光谱波段数目多等特点，逐渐成为领域内的研究热点(张良培和李家艺，2016)。但是由于成像光谱仪空间分辨率的限制和地物的复杂多样性，一个像元中通常包含多种地物，该像元称为混合像元，混合像元普遍存在于高光谱遥感数据中(陈晋 等，2016)。混合像元现象制约了高光谱遥感地物模式识别和解译等高层应用。因此，如何分解高光谱混合像元(高光谱解混)成为高光谱数据处理领域的一个关键问题。其目的是将混合像元分解为一组纯物质(端元)和对应的成分比例(丰度)(Bioucas-Dias等，2012)。

线性混合模型(linear mixing model, LMM)假设端元之间不存在相互作用，混合像元由一组端元和对应成分的丰度线性组合而成，是常见光谱解混算法的物理假设，简单高效且物理意义明确。基于LMM的高光谱解混方法包括基于几何学的方法、基于统计学的方法和基于稀疏回归的方法(朱玲 等，2020)。

基于几何学的方法认为，高光谱数据中的每一个像元都对应多维特征空间中的一点，所有像元都位于多维特征空间内的一个单纯形几何体中，而单纯形的顶点为各端元，因此可以通过求解单纯形的顶点得到端元光谱(蓝金辉 等，2018)。纯端元假设指高光谱数据中至少存在1个像元是只由1种物质组成的。根据是否采用纯端元假设，基于几何学的方法大致可以分为两类。一类是假设高光谱数据中存在纯端元，代表性工作有N-FINDR(Winter，1999)，顶点成分分析(vertex component analysis，VCA)(Nascimento和Dias，2005)；另一类是最小体积法，适用于没有纯端元的数据，代表方法有最小体积单纯形分析(minimum volume simplex analysis，MVSA)(Li和 Bioucas-Dias，2009)和采用分裂与拉格朗日增广的单纯形辨识方法(simplex identification via splitting and augmented Lagrangian，SISAL)(Bioucas-Dias，2009),这类算法试图寻找包含所有像元的最小体积单纯形来提取端元。

上述方法需要先提取端元，再反演丰度。基于统计学的方法可以同时获取端元和丰度,如非负矩阵分解(nonnegative matrix factorization，NMF)方法(宋晓瑞 等，2020)。由于NMF的解不唯一，为了克服这一特性，需要引入约束，如丰度重加权稀疏NMF(reweighted sparse NMF，ARSNMF)(贾麒 等，2020)引入了重加权1约束，稀疏和正交约束NMF(sparse and orthogonal constrained NMF，SONMF)(陈善学和储成泉，2019)引入了1/2约束和正交性约束。

随着美国地质调查局(United States Geological Survey，USGS)光谱库的出现，基于稀疏回归的方法受到业界关注(袁静 等，2018)。该类方法将事先收集的大量纯物质光谱组成光谱库，并作为端元字典，避免端元提取不准确带来的影响，进而通过稀疏回归模型进行求解。Bioucas-Dias和Figueiredo(2010)对丰度矩阵施加非负性和稀疏性约束，提出了基于变量分裂和增广拉格朗日的稀疏解混模型(sparse unmixing by variable splitting and augmented Lagrangian，SUnSAL)，利用拉格朗日方法和交替方向乘子法对模型进行求解。由于光谱库中的端元数量远大于高光谱数据中的端元数量，因此丰度矩阵中会有少量的非零行，体现出行稀疏性。利用丰度矩阵的行稀疏性，Iordache等人(2014)在SUnSAL模型的基础上，将丰度矩阵的1范数换成2,1范数，提出了协同稀疏回归模型(collaborative SUnSAL，CLSUnSAL)。在此基础上，Shi等人(2018)提出使用2,0范数刻画丰度矩阵的行稀疏性，在保持高精确度的情况下得到了更稀疏的解。高光谱图像中邻近像元对于同一种端元可能具有相似的丰度，因此Iordache等人(2012)用全变分(total variation，TV)刻画像元的局部同质性和分段平滑性，提出了基于全变差和变量分裂增广拉格朗日的稀疏解混模型(sparse unmixing via variable splitting augmented Lagrangian and total variation，SUnSAL-TV)。为了进一步挖掘空间信息，黄伟等人(2020)提出了基于局部加权低秩先验的高光谱稀疏表示解混方法(hyperspectral sparse unmixing based on local weighted low-rank prior，HSULWLrP)，使用加权低秩先验挖掘局部块的低维结构特征。上述模型都假设不同波段高光谱数据的噪声服从高斯分布且噪声强度相同。然而,实际高光谱数据噪声复杂，可能包含脉冲、死线等其他噪声。基于这一认识，Li等人(2020)提出了分波段的稀疏解混方法(sparse unmixing method with the bandwise model，SUBM)。

本文对高光谱图像进行非负稀疏分量分解建模，提出了一种基于非负稀疏分量分析的鲁棒解混方法(robust unmixing for hyperspectral images with sparse component analysis，RUnSCA)。RUnSCA在最大后验框架下，考虑高光谱数据不同波段的混合噪声特性，建模为稀疏性结构噪声，并根据其稀疏性施加1范数约束，同时根据丰度的行稀疏性为其施加2,0范数约束，最后加入TV正则项以促进丰度图的平滑性。在模拟数据集和真实数据集上的实验均表明，RUnSCA对混合噪声鲁棒，可以有效克服同批光谱数据之间的波形形态结构差异，在混合噪声条件下，RUnSCA可以获得最优的信号与重构误差比(signal to reconstruction error，SRE)。

1 相关工作

1.1 线性混合模型

线性混合模型(LMM)(张兵，2016)认为混合像元的光谱是若干端元的光谱按照相应的比例系数线性混合得到的，其数学描述为

(1)

式中，y∈RB表示高光谱数据中单个混合像元的光谱向量，B为波段数量，ai∈RB表示端元光谱向量，xi表示ai在混合光谱y中所占的比例系数，即丰度，M表示y中包含的端元数量，n表示噪声。LMM的矩阵形式可以表示为

Y=AX+N

(2)

式中，Y∈RB×P表示高光谱数据的矩阵形式，包含B个波段，P个像元，其中每一列为单个像元的光谱向量，A∈RB×M表示端元字典，包含M个端元，X∈RM×P表示丰度，N∈RB×P表示噪声。

一般地，丰度通常满足“非负性”约束以及“和为一”约束(Heinz和Chang，2001)，即对于每个端元的丰度xi，满足

(3)

1.2 稀疏解混

稀疏解混将由大量纯端元光谱组成的光谱库作为端元字典，然后从光谱库中挑选合适的端元，再计算对应的丰度(Iordache等，2011)。因为光谱库中的端元数量远大于高光谱数据中的端元数量，使得丰度系数矩阵是稀疏的，因此可以针对其稀疏性施加约束，得到约束的0优化问题，具体为

(4)

(5)

2 本文方法

2.1 非负稀疏分量分解建模

大多数稀疏解混方法都假设高光谱数据所有波段的高斯噪声强度相同，然而在实际情况中，不同波段的高斯噪声的强度不同，并且通常被其他类型的噪声污染(Zhang等，2014)，如脉冲噪声、死线噪声等，这类噪声可以归类为稀疏噪声。在混合噪声的情况下，对高光谱数据的第i个波段建模，具体为

Yi=(AX)i+Si+Ni,i=1,…,B

(6)

综上，数据似然概率可表达为

(7)

式中，c为常数。Wi,i=1/σi(i=1,…,B)∈RB×B可以视为一个权重矩阵，高斯噪声的方差越大，波段的权重越小，从而考虑到高斯噪声强度因波段而异的情况。然后在最大后验框架下，估计丰度和稀疏噪声，建立非负稀疏分量分解优化模型，具体为

(8)

式中，第2个等式根据贝叶斯定理得到，对第2个等式应用-ln(·)后可以得到第3个等式，由于c为常数，在优化问题中可以将lnc这一项舍弃，得到最后一个等式，p(AX+S)表示AX+S的先验分布，可以看做施加在丰度和稀疏噪声上的先验知识约束。

由于光谱库中端元的数量远大于高光谱数据中端元的数量，因此丰度矩阵中会有少量的非零行，体现出全局行稀疏性(Iordache等，2014)。利用丰度矩阵的全局行稀疏性，可以为丰度矩阵施加2,1混合范数约束，从而更加有效地刻画丰度矩阵的稀疏性。此外，高光谱数据中被稀疏噪声S污染的像元较少，体现出稀疏性，可以用0范数对其建模。由于0范数的优化问题NP难，因此对S施加1范数约束。进而，式(8)可以表示为

(9)

(10)

由于光谱的可变性(Bateson等，2000)，丰度的“和为一”约束在实际情况中一般难以满足，所以本文只考虑“非负性”约束。高光谱图像中，局部相邻的像元通常包含相似的物质，其光谱也具有相似性，在图像中表现出来就是分段平滑性，由于丰度矩阵每列中的元素代表的是相应物质在像元中所占的成分比例，因此可以认为丰度矩阵也是分段平滑的(Iordache等，2012)。式(10)加入TV正则项，得到基于非负稀疏分量分析的高光谱鲁棒解混模型，具体为

(11)

式中，λTV>0表示正则化参数，TV(X)表示各向异性TV(anisotropic TV)。

令Hh:RM×P→RM×P表示丰度X相邻像元之间水平差分的线性算子(Iordache等，2012)，即HhX=[d1,d2,…,dP]，其中di=xi-xih，xih为xi的水平邻像元。令Hv:RM×P→RM×P表示丰度X相邻像元之间垂直差分的线性算子，即HvX=[v1,v2,…,vP]，其中vi=xi-xiv，xiv为xi的垂直邻像元。在上述两种算子的基础上，本文定义

(12)

因此模型(11)可以写为

(13)

2.2 模型求解

引入辅助变量V1,V2,V3,V4,V5后，模型式(13)可以写为

(14)

s.t.V1=S;V2=X;V3=X;V4=HV3;V5=X

式中，lR+(V5)是一个指示函数(indicator function)，当V5≥0时，lR+(V5)=0，否则为+∞。需要注意的是，式(13)中V4=HV3与其他约束不对称，这样可以将两个变量的优化解耦。式(14)可以写成紧凑形式，具体为

(15)

式中

(16)

式(15)的增广拉格朗日方程为

(17)

式中，D/μ为拉格朗日乘子，D=(D1,D2,D3,D4,D5)。

上述约束最小化问题是多变量联合最小化问题，直接联合最优复杂度较高。本文采用交替方向乘子法(alternating direction multiplier method，ADMM)将其分解为若干个子问题，对每一个子问题交替迭代求解。首先固定其他变量，更新变量X，具体为

(18)

对式(18)求导，可以得到X的解，即

(19)

(20)

其解可以用软阈值法(Bioucas-Dias，2009)求得，即

(21)

式中，Sτ(δ)=sgn(δ)max(|δ|-τ,0)表示软阈值收缩算子，sgn(δ)表示指示函数，τ表示阈值。

更新V1时的优化问题为

(22)

其解为

(23)

更新V2时的优化问题为

(24)

其解可由行硬阈值法(Shi等，2018)求得，即

(25)

式中，RHτ(Φ)表示行硬阈值(row hard threshold)函数，定义为

(26)

式中，τ>0为阈值，Φ(i,:)表示矩阵Φ的第i行。

更新V3时的优化问题为

(27)

其解为

(28)

式(28)可以通过快速傅里叶变换(fast Fourier transform, FFT)(Sun等，2018)求解，具体为

(29)

更新V4时的优化问题为

(30)

其解也可以用软阈值法求得，即

(31)

更新V5时的优化问题为

(32)

其解为

(33)

迭代过程中，本文采用基于最小误差的高光谱信号识别(hyperspectral signal identification by minimum error，HySime)(Bioucas-Dias和Nascimento，2008)估计高斯噪声，进而计算得到W。

用ADMM求解RUnSCA的算法流程如下：

输出：X,S。

2) 估计Y-S(k)的高斯噪声，更新W；

3) 通过式(19)更新X(k+1)；

4) 通过式(21)更新S(k+1)；

15)k=k+1；

16) end；

17) 返回X=X(k+1),S=S(k+1)。

在上述算法流程中，计算较复杂的步骤是W、X和V3的求解，它们的复杂度分别为O(B2P)、O(B2P)和O(BPlogP)，其中B为光谱波段数量，P为像元数量。鉴于B通常大于logP，因此RUnSCA的总体计算复杂度为O(B2P)。

3 实验结果与分析

3.1 数据集

为测试方法的有效性，采用3组数据集进行实验。

2)模拟数据集2(data cube 2, DC2)。该数据集的生成参考了Miao和Qi(2007)的方法，包含224个波段，共64×64个像元，不含纯像元。

3.2 实验设置

实验将RUnSCA与分波段稀疏解混方法SUBM(Li等，2020)、基于TV正则化的方法SUnSAL-TV(Iordache等，2012)、RSSUn-TV(Wang等，2019)、基于协同稀疏的方法CLSUnSAL(Iordache等，2014)和CSUnL0(sparse-cooridinated hyper-spectral un-mixing using0norm)(Shi等，2018)等5种具有代表性的方法进行对比。

采用信号与重构误差比(SRE)以及均方根误差(root mean square error， RMSE)作为性能评价指标，具体定义为

(34)

(35)

实验参数固定μ=10-2，ε=10-6，N=103，正则化参数λ,α和λTV从集合中取值，即λ,α,λTV∈{10-5,10-4,…,10-1,100,101,…,10-4,105}。实验时考虑这些参数的所有组合，采用最优的性能评价指标值作为方法的最终性能。

实验环境为MATLAB R2015b软件，操作系统为64位Windows10，硬件平台配置为Intel Core i7-9700K CPU，主频为3.60 GHz，内存为16 GB。

3.3 模拟数据集1(DC1)实验

在DC1数据集上进行光谱解混实验。不同方法在DC1数据集上解混得到的丰度图如图1所示。可以看出，CLSUnSAL、CSUnL0、SUnSAL-TV和RSSUn-TV的丰度图中有明显的噪声。SUBM和RUnSCA两列背景最为纯净，二者相较，RUnSCA的丰度图噪声更少，与真实丰度图更为接近。

图1 不同方法在DC1数据集上得到的丰度图

3.3.1 超参数分析

为了确保原始残差与对偶残差之间的比率维持在一定区间内，初始参数值μ=10-2，迭代过程中不断更新参数μ。图2为RUnSCA在DC1数据集上使用不同正则化参数λ、α和λTV的解混性能结果。图2(a)为RUnSCA在不同的正则化参数λ下的SRE。可以看到，当λ<101时，λTV=10-2的SRE明显优于λTV=102的SRE，说明较大的λTV会降低RUnSCA的解混精度；当λ>101时，4条曲线接近一致。设置μ=10-2,λ=100，不同的α和λTV对RUnSCA的影响如图2(b)(c)所示。可以看出，较小的λTV带来更高的SRE，在α>10-1或λTV>101时，SRE保持不变。

图2 RUnSCA在DC1数据集上使用不同的正则化参数λ、α和λTV得到的SRE

3.3.2 多种噪声强度下的有效性

对DC1数据集加噪，改变高斯噪声的强度，可得到平均信噪比不同的数据集。图3为不同信噪比条件下各种方法对比。

图3 不同方法在不同信噪比条件下在DC1数据集上得到的SRE

从图3可以看到，随着信噪比的增大，所有方法的SRE都随之提高，但RUnSCA的SRE始终高于其余对比方法。

3.4 模拟数据集2(DC2)实验

在DC2数据集上进行实验。图4为混合噪声条件下不同方法在DC2数据集上解混的丰度图。图4中第1列为5种端元的真实丰度图。可以看出，CLSUnSAL、CSUnL0、SUnSAL-TV和RSSUn-TV的丰度图与真实丰度图差距较大，SUBM和RUnSCA的结果与真实丰度图十分接近，说明考虑各波段高斯噪声的不同强度以及稀疏噪声建模可以有效提升解混性能。总的来说，与其他方法相比，RUnSCA的结果最接近真实丰度图。

表1展示了RUnSCA在DC2数据集上使用不同参数得到的SRE。可以看到，在λ=0和λTV=0的情况下，RUnSCA获得了更低的SRE，证明了2,0范数和TV正则项的有效性。

表1 RUnSCA在DC2数据集上使用不同参数得到的SRE

1)利用丰度矩阵的低秩性(Zhao和Yang，2015)，本文将式(11)中的2,0正则项替换为低秩正则项，使用加权核范数求解相关子问题，该方法称为RUnSCA-Rank。

2)由于稀疏噪声存在于高光谱数据中的少数波段，因此可以认为稀疏噪声也具有行稀疏性。本文将式(11)中的1正则项替换为2,0正则项，该方法称为RUnSCA-Row。

3)各向同性TV(isotropic TV)与本文采用的各向异性TV(anisotropic TV)分别计算为

(36)

(37)

本文将式(11)中的各向异性TV替换为各向同性TV，该方法称为RUnSCA-Iso。

图5展示了上述3种方法与RUnSCA在混合噪声条件下的对比。可以看到，RUnSCA-Iso与RUnSCA的结果十分接近，在高信噪比情况下几乎一致。但从整体上看，RUnSCA-Iso仍然微弱于RUnSCA，说明在混合噪声情况下，各向异性TV优于各向同性TV，但差距不大。此外，在低信噪比情况下，RUnSCA-Rank的SRE接近RUnSCA，RUnSCA-Row的SRE明显弱于RUnSCA，然而在高信噪比时情况相反，说明低秩正则项能有效抑制噪声，但无法取得比2,0正则项更稀疏的结果，并且2,0正则项对高强度噪声较为敏感。总体来看，RUnSCA方法优于其余对比方法。

图5 不同方法在不同信噪比条件下在DC2数据集上得到的SRE

3.4.3 不同类型噪声下的有效性

表2为不同方法在不同噪声情况下在DC2数据集上得到的SRE。可以看出，当高光谱数据仅被高斯噪声污染时，RUnSCA获得了最高的SRE，这是因为其考虑了高斯噪声在不同波段上的变化。当只被脉冲噪声或死线噪声污染时，RUnSCA的SRE同样最高，原因在于RUnSCA考虑了脉冲、死线等稀疏噪声，同时针对它们的潜在稀疏性进行建模，证明了对稀疏噪声建模可以有效提高解混的精度。可以看到，基于TV正则化的方法(SUnSAL-TV、RSSUn-TV)获得的SRE比基于协同稀疏的方法(CLSUnSAL、CSUnL0)高，说明TV正则化具有更好的抗噪性能。此外，RSSUn-TV获得了比SUnSAL-TV更高的SRE，因为除了采用TV正则化之外，RSSUn-TV还采用2,0范数来表示丰度的全局行稀疏性。然而，在大多数情况下，CSUnL0的SRE比CLSUnSAL低，这说明与2,1范数相比，2,0范数对高强度的噪声更加敏感。在同时存在高斯、脉冲和死线噪声时，RUnSCA的SRE最高，说明RUnSCA的抗噪性能最好。总结来看，在不同条件下，RUnSCA都能得到最高的SRE，说明其对混合噪声具有鲁棒性。

表2 不同方法在不同噪声情况下在DC2数据集上得到的SRE

3.4.4 收敛性分析

图6为不同方法在DC2数据集上的收敛曲线。可以看到，所有方法都在100次迭代以内快速收敛，然后趋于稳定。尽管RUnSCA的收敛曲线有小幅波动，但整体来讲，收敛效果比其余方法更好。

图6 不同方法在DC2数据集上的收敛曲线

3.4.5 单条光谱解混的有效性

提取DC2数据集中单条光谱像元曲线作为输入数据，即y∈R224。图7为原始光谱、含噪光谱以及不同方法重构后的光谱，表3为不同方法重构该光谱得到的RMSE。可以看出，图7(h)中RUnSCA重构后的光谱最接近原始光谱，表3中RUnSCA的RMSE最小，综合说明RUnSCA能较好地去除高斯和稀疏噪声，重构性能较好。

表3 不同方法重构单条光谱得到的RMSE

图7 DC2数据集中单个像素的原始光谱、含噪光谱以及不同方法的重构结果

3.5 真实数据集实验

图8显示了不同方法在Cuprite数据集上得到的3种矿物alunite、buddingtonite、chalcedony的丰度图。图8(a)是用Tetracorder软件得到的3种矿物的分类图，其将每个像素分类为特定矿物，这里仅用做比较分析。图8(g)是RUnSCA的解混结果，色调越热表示该矿物的含量越高。可以看出，SUBM和RUnSCA得到的丰度图相比其他方法蕴含更多细节。从视觉效果上看，RUnSCA得到的丰度图具有可比性，证明RUnSCA是有效的。

图8 不同方法在Cuprite数据集上得到的3种矿物的丰度图

图9显示了RUnSCA对Cuprite数据集上单条混合像元光谱曲线进行解混和重构的结果。图9(a)为该混合像元的光谱以及使用RUnSCA重构的光谱。图9(b)为该混合像元解混后的前5个端元，分别为Cheatgrass、Hydroxyl-Apatite、Kaolin/Smect、Alunite和Erionite+Merlinoit，对应的丰度系数分别为0.131 0, 0.091 0, 0.073 2, 0.053 2和0.052 9。

图9 对Cuprite数据集上单条像元光谱解混和重构的结果

不同方法在DC1、DC2和Cuprite数据集上迭代100次的运行时间如表4所示。由于需要进行差分计算，因此SUnSAL-TV、RSSUn-TV和RUnSCA的时间最长，其中RUnSCA的运行时间最多。总体而言，运行时间在可接受范围内。

表4 不同方法在DC1、DC2和Cuprite数据集上的运行时间

4 结 论

考虑到高光谱数据中存在混合噪声，且不同波段高斯噪声强度不同，本文对高光谱图像进行非负稀疏分量分解建模。该模型基于最大后验概率框架，对稀疏结构性噪声建模，同时基于丰度的行稀疏性与相邻像素之间的局部空间平滑特性，用2,0促进丰度的行稀疏性，TV正则项促进丰度图的分段平滑。最后采用ADMM优化算法迭代求解，设计了一种新的解混方法。在模拟数据集和真实数据集开展了收敛性、鲁棒性和单条光谱解混等多组实验。实验结果表明，与代表性解混方法相比，提出方法在混合噪声条件下表现出了优异的解混效果。同时，实验综合分析了丰度的行稀疏性与全局低秩性、各向异性TV与各向同性TV，以及稀疏噪声的1范数与2,0范数对解混性能的影响，表明本文模型具有较好的噪声稳健性和盲源分离性能。应该看到，本文研究是聚焦于在线性混合模型假设下的结果，因此一方面可以进一步考虑拟线性和非线性混合假设下的理论与算法推广；另一方面可以有效结合数学模型解混和深度先验学习，进一步挖掘高光谱图像的非线性复杂结构特征，从而设计模型启发的高光谱解混深度网络算法。