基于L0稀疏约束的近似NMF高光谱解混

刘雪松，谭文群，段卓镭，张志鹏，彭天亮

(南昌工程学院 江西省水信息协同感知与智能处理重点实验室，江西 南昌 330099)

高光谱遥感图像具有很高的光谱分辨率，能够在很窄的光谱波段上获得连续的光谱影像信息，因而被广泛用于环境监测、农业生产、食品检测以及地质勘探等领域[1]。高光谱遥感影像所获取的像元一般是混合像元，因此需要利用高光谱解混技术对混合像元进行分解，从中提取出相应的地物信息(端元)，及这些端元所占像元的比例(丰度)。如何精确而又方便的处理混合像元分解问题，这是高光谱遥感技术所要处理的关键问题。

非负矩阵分解(NMF)旨在将非负矩阵X分解为非负矩阵A和S的乘积，即X≈AS，这种模型能够很好的描述高光谱解混的线性模型。高光谱数据解混通常还要加入两个限制项：端元丰度非负性约束(ANC)和丰度和为一约束(ASC)，但这两个限制项只是先验知识的一部分。由于NMF经典理论是基于整块数据进行分解，当数据很大时，计算量也非常大，而高光谱数据往往数据维数特别大[2]。并且，光谱库中的反射光谱曲线个数远远大于实际有效的端元个数，这就使得观测的每一个混合像元反射光谱曲线在光谱库下的丰度系数向量可以被看成是“稀疏”的[3]。这就为“稀疏表示”的思想引入到高光谱混合像元分解中奠定了基础。

近年来，基于稀疏性约束的高光谱解混是研究的热点。该研究思路主要分为两类，一类是基于光谱库的稀疏性高光谱解混，另一类是基于非负矩阵分解的稀疏性高光谱解混。本文主要讨论的是基于非负矩阵分解的稀疏性高光谱解混，其中代表性的工作有：贾森[4]等提出添加分段光滑约束和稀疏度约束的非负矩阵分解用于高光谱解混，稀疏性约束中最典型的是L0范数稀疏，其稀疏性效果很好，但因存在NP难题，求解复杂，所以实际应用很难。Candes[5]等提出利用L1范数替代L0来进行求解，这样使得算法的速度加快，但是稀疏性却减弱了。钱沄涛等利用L1/2范数约束端元对应丰度系数的非负矩阵分解算法(L1/2NMF)进行解混[6]。

为了保证稀疏特性的稳定和算法的快速性，本文在L0正则化理论的基础上对基于NMF的高光谱解混方法进行了改进，提出了一种近似L0NMF框架(NMFL0_A，NMFL0_S)的算法应用于高光谱解混，通过对端元矩阵和丰度矩阵的联合稀疏约束来提升NMF算法的稀疏表达能力。相比于经典的NMF算法，近似NMF可以利用部分数据集解混出相关的端元和丰度的信息，并且解混的精度和运算速度明显提升。

1 稀疏正则化NMF

1.1 L0稀疏正则化NMF光谱解混技术

在经典NMF模型中，由于其对初始化数据非常敏感。所以有学者提出了许多关于稀疏性正则化的NMF理论，理论上，0

(1)

式中Sn(k)是图像n像素处的k端元的丰度分数。在稀疏正则化的NMF理论中，最典型的就是L0范数稀疏。其目标函数描述如下：

(2)

式中X∈RM×N为高光谱图像的数据矩阵；A∈RM×K为包含了K个端元的端元矩阵；S∈Rk×N为丰度系数矩阵；α‖S‖0是 L0稀疏的附加约束函数；α∈R为正则化项，用来表示稀疏特性的附加约束函数的权重。

尽管L0范数的稀疏性很好，但是L0正则化器在高光谱解混中是一个NP问题，在实际应用中无法解决。因此有学者对L0正则项提出相关逼近算法，Tao和Candes证明了在RIP条件下可以用L1范数代替L0范数求解，其原理是在RIP条件下，L1范数和L0范数具有相同的解，且L1范数比L0范数求解简单，速度更快[8]。所以L0范数在稀疏解混上远远没有用L1范数求解方便。

1.2 L1/2稀疏正则化NMF光谱解混技术

虽然高光谱解混中的稀疏特性和L1正则化器都已经得到了广泛的使用，但是为了得到稀疏效果更好的解混方法，Qian和Jia将L1/2范数应用到了非负矩阵分解算法中来[9]，即L1/2NMF。它在光谱解混中的计算更加简单，并且产生比基于L1范数的算法更稀疏，更准确的结果。因为，L1/2情况与式(1)～(2)情况类似，所以归纳之后，得到基于L1/2NMF的目标函数描述如下：

(3)

(4)

对(3)求偏导，得

(5)

(6)

根据梯度下降法得

A←A-ξA(ASST-XST).

(7)

(8)

式中ξA和ξS为迭代步长。为了保证非负性，可将式(7)～(8)式变形为

A←A[(XST)/(ASST)].

(9)

(10)

当迭代A和S达到规定的停止条件时，更新迭代停止。基于L1 /2稀疏性非负矩阵分解的高光谱混合像元分解算法步骤如下：

①使用HySime算法来估算高光谱中端元个数p；

②根据X上的稀疏度，测量估算稀疏的权重因子α；

③利用VCA-FCLS对A和S进行初始化，并且对S归一化处理；

④重复利用式(9)～(10)，对矩阵A和S每行元素进行更新，直到达到最大迭代次数。

2 基于L0稀疏约束下的近似NMF

目前，基于稀疏性高光谱混合像元分解技术是利用已知的光谱库作为端元集合来进行分解[1]。这一过程中，通常是对丰度矩阵S加入稀疏性约束，以此作为附加条件来进行高光谱混合像元分解。虽然这种方法解混出来的效果也比较理想，但是，此方法解混出来的端元矩阵A却不是稀疏矩阵。并且这种稀疏约束的方法在高光谱解混中比较单一。所以，为了让矩阵A和S同时具有稀疏约束的效果，并且解混的精度更好。在L0正则化理论的基础上对基于NMF的高光谱解混方法进行了改进，并与无约束的NMF技术(乘法更新规则、交替非负最小二乘)结合起来，提出了一种近似L0NMF框架(NMFL0_A，NMFL0_S)，即联合稀疏约束模型，并将其应用到高光谱遥感图像混合像元分解中，从而实现了对端元矩阵和丰度矩阵的联合稀疏约束。

2.1 通过乘法更新的NMF

NMF算法通常以两个阶段的迭代方式进行，即交替更新A和S，从而让非凸优化问题变成凸优化问题，方便进行求解。这也是本文所要关注的矩阵更新规则。Lee and Seung[10]展示了通过乘法更新的NMF。乘法更新规则如下：

S←S⊗[(ATX)/(ATAS)].

(11)

A←A⊗[(XST)/(ASST)].

(12)

式中“⊗”和“/”分别表示按元素乘法和除法。显然，因为X的元素是非负值的，所以这些更新规则保持了A和S的非负性。

此外Paatero and Tapper[11]也进行过NMF算法的更新，他们通过交替的非负最小二乘进行A和S的交替更新，更新规则如下

(13)

(14)

这里注意，‖X-AS‖F在A或S中分别是凸的，而在A和S整体中则是非凸的。通过求解具有多个右手边(即X的每一列一个)的非负约束最小二乘问题(NNLS)，可以找到每个子问题的最佳解决方案。因此，此方案也称为交替非负最小二乘(ANLS)。

2.2 具有L0稀疏约束的NMF

目前，对于L0稀疏约束的NMF方法还是相对较少的。其中，K-SVD算法[12]可以被认为是NMF算法在S列上具有L0稀疏约束；另外，在文献[13]中，也提出了一种在S列上具有L0稀疏约束的NMF算法，在稀疏编码阶段，他们使用了一个非负的最小角度回归和选择LARS[14]，简称为NLARS，也可以被认为是具有L0稀疏约束的NMF算法。现在，要介绍的是在A和S列上都具有L0稀疏约束的NMF的求解方法，形式上，由下面式子表示：

(15)

(16)

式(15)～(16)分别称为NMFL0_S和NMFL0_A，参数L是Si和Ai中非零条目的最大允许个数。对于NMFL0_S的稀疏约束意味着X中的每一列都用一个极大值的锥形组合来表示L个非负基向量；而NMFL0_A，稀疏约束在有限的支持下强制执行基向量，比如，如果X的列包含图像数据，那么A上的稀疏约束鼓励采用基于部分的表示。

2.3 NMFL0_S

由式(11)～(14)可知，NMF算法以两阶段的迭代方式进行，所以在L0稀疏约束的NMF算法中，对NMFL0_A和NMFL0_S也应用了同样的原理，NMFL0_S算法过程由算法1给出。NMFL0_S算法步骤如下：

①随机初始化基础矩阵A；

②设置迭代次数i；

③非负稀疏编码：固定基础矩阵A，对X进行非负稀疏编码，固定基础矩阵A，导致稀疏，得到稀疏的非负矩阵S；

④矩阵更新：保持S的稀疏结构，更新A和S。

在针对NMFL0_S时，它的第一阶段是先进行非负稀疏编码，但稀疏编码问题是NP难的，所以需要一个近似值。为此，Robert等将正交匹配追踪(OMP)和非负最小二乘(NNLS)结合，提出了稀疏非负最小二乘(sNNLS)[15]。从而将凸问题转化成了优化的非凸问题，得到每个子问题的最优解。为了得到更好的稀疏效果，进一步提升算法的计算速度，在进行匹配追踪时，提出了一种新的匹配追踪原则，将非负最小二乘与反向匹配追踪结合，即反向匹配追踪(rsNNLS)。它不仅适用于非负的框架，而且适用于匹配追踪原则。第二阶段是矩阵更新阶段，这一阶段的目标是增强基础矩阵A。在这一步骤中，可以适当调整S中的非零值，但是不允许零值变成非零值，即要保持S的稀疏结构。

NNLS算法步骤如下：

①由Z={1，…，k}，p为非零数，令s=0，r=x-As，α=ATr；

②当|Z|>0和∃i∈Z：αi>0时，有i*=arg maxα；

③根据Z←Zi*，P←P∪i*，P←{i|si>0}；

⑤当∃j∈p：zj<0时，有α=minsk/(sk-zk)其中，k∈p；

⑥根据s←s+α(Z-s)，Z←{i|si=0}，P←{i|si>0}；

⑧结束片刻；

⑨根据s←Z；

⑩得到r=x-As，α=ATr。

算法中x和s分别是数据矩阵和系数矩阵的列向量，Zp代表只含有集合P中元素的矩阵。活动集Z和非活动集P不相交的包含h个条目的索引。具有Z中索引的条目保持为零，其非负性约束是可变的，而含有P索引的条目则保持正值。在此算法的外部循环中索引从Z移动到P，直至获得最优解[16]。步骤⑤到⑦是内部循环，纠正了向量中存在的负解。

反向稀疏(rsNNLS)算法步骤如下：

①由s=NNLS(x，A)，得到Z={i|si=0}，P={i|si>0}；

②当‖s‖0>L时，有j=arg minsi，其中ip；

③根据sj←0，Z←Z∪{j}，P←P{j}；

④重复NNLS算法中的④～⑨步骤；

上述算法是从步骤①中的最佳非稀疏解开始的，虽然L0范数的解大于L，但是解向量中的最小条目被设置为零，并且其索引是从无效集合P移到有效集合Z(步骤②～③)。随后，使用NNLS算法中的④～⑨步骤，在NNLS的意义上用P中的剩余基向量对数据向量进行近似。

2.4 NMFL0_A

算法框架设计的想法来源于Hoyer的L1NMF算法[14]。在此，依旧遵循前面矩阵更新的两个阶段的迭代方法。首先固定系数矩阵S，然后计算基础矩阵A的最优无约束解。接下来，将基础矩阵向量投影到欧式空间中最接近的非负向量上，以满足L0约束的条件。这一步很容易，只需保留向量中L个最大项，其余项置零。在矩阵更新阶段，只需调整A中的非零条目，即让非零项变为零，从而维持A的稀疏结构[17]。NMFL0_A算法步骤如下：

①随机初始化系数矩阵S；

②设置迭代次数i；

③得到AT=NNLS(XT，ST)；

④令j=1：k；

⑤将Ai中的D-L最小值设置为零；

⑥矩阵更新：增强A和S，并保持A的稀疏结构。

3 实验结果与讨论

实验选取真实高光谱数据库，对经典的高光谱解混及其约束下的算法(NMF，L1/2NMF)及基于改进L0稀疏约束下的近似NMF算法(NMFL0_A，NMFL0_S)光谱解混性能分别进行测试和比较，从光谱角距离(SAD)、均方根误差(RMSE)、和运算速度这三方面来评价解混的精度和效果。

(17)

(18)

3.1 仿真数据实验

第一个仿真数据，采用高光谱数据库Samson。Samson是一个简单的数据集。在这个图像中，有952×952像素。每个像素记录在156个波段中，覆盖波长范围是401～889 nm。光谱分辨率最高可达3.13 nm。由于原始图像太大，计算成本非常昂贵，使用了95×95像素的区域。它从原始图像中的(252，332)像素开始。这些数据不会被空白信道或糟糕的信道所破坏。具体来说，这幅图中有3个目标：“#1 Rock”“#2 Tree”和“#3 Water”。

本实验是在Samson数据集上，分别执行了NMFL0_S和NMFL0_A两种算法。实验所得数据均进行了100次计算，取平均值。端元数P为3，k为3，在NMFL0_S和NMFL0_A中L分别为3和156，参数i分别为5和30。表1为相应算法下Samson高光谱数据进行混合像元分解的SAD和RMSE的精度值。从表1中可以得出：NMFL0_A算法相对其他算法SAD值降低0.01～0.04，NMFL0_S算法相对其他算法RMSE值降低0.11～0.12，解混的精度明显提高。并且从表1中可以看出两种近似L0稀疏约束算法与其他两种算法相比，NMFL0_S在运算速度上也有所提升，优势明显，NMFL0_A运行速度稍慢。

表1 4种算法对Samson高光谱数据进行混合像元分解的精度和速度的均值比较(粗体数字表示最好结果)

为了更清晰的比较几种算法解混的精度，本文还将几种算法的端元光谱图像和丰度图像进行了比较分析，如图1～2所示。图1中给出了几种算法(NMF，L1/2NMF，NMFL0_A)对Samson高光谱数据进行的混合像元分解后的端元估计值与真实端元值的比较图，从图(1)中可以看出，三种算法的端元光谱基本都接近真实的光谱，相差不大。在Tree和Water的光谱曲线中，NMFL0_A的鲁棒性最好，优势非常明显，最接近真实光谱曲线，拟合效果最好。其次是L1/2NMF算法，最差的是NMF，在Soil中，偏离真实光谱曲线较大。

图1 基于几种算法下的端元信息图

图2是几种算法(NMF，L1/2NMF，NMFL0_S)对Samson高光谱数据进行的混合像元分解后的端元对应的丰度图像。从图2中可以看出三种算法都能较好的接近真实的丰度图。在树和水部分，NMFL0_S算法的丰度分布图明显与真实的丰度图更加接近，而且图像在水部分，图(l)相比于L1/2NMF算法的丰度图(k)要清晰很多，而NMF算法，则存在像元丰度相互交融的情况。在树部分图(h)相比于NMF算法的丰度图(f)要更接近真实端元的丰度信息。所以综合图1～2，可以看出NMFL0_S和NMFL0_A这两种算法相比NMF和L1/2NMF算法更加优越，效果更好。

3.2 仿真数据实验

第2个仿真数据，采用Jasper高光谱数据库，该数据库包含512×614像素，每个像素记录在范围从380～250 0nm的224个通道中，光谱分辨率高达9.46nm。这个高光谱图像太复杂，因此，考虑100×100像素的子图像。在移除通道1～3，108～112，154～166和220～224后(由于密集的水蒸气和大气效应)，保留了198个通道(这是光谱解混分析的常见预处理)。在这些数据中有4个目标：“#1 Road”，“#2 Soil”，“#3 Water”和“#4 Tree”。

本实验是在Jasper高光谱数据库上，分别执行了NMFL0_S和NMFL0_A两种算法。实验所得数据均进行了100次计算，取平均值。端元数P为4，参数i为5，k为4，在NMFL0_S和NMFL0_A中L分别为4和198。表3为相应算法下Jasper高光谱数据进行混合像元分解的SAD和RMSE的精度值。从表2中可以得出：NMFL0_A算法相对其他算法SAD值降低0.04～0.11，NMFL0_S算法相对其他算法RMSE值降低0.05～0.08，解混的精度明显提高。并且从表2中可以看出两种近似L0稀疏约束算法与其他两种算法相比，在运算速度上有了大幅度提升，优势明显。

表2 4种算法对Jasper高光谱数据进行混合像元分解的精度和速度的均值比较(粗体数字表示最好结果)

(a)真实(沙岩)丰度图 (b)NMF(沙岩)丰度图 (c)L1/2NMF(沙岩)丰度图(d)NMFL0_S(沙岩)丰度图 (e)真实(树)丰度图 (f)NMF(树)丰度图 (g)L1/2NMF(树)丰度图 (h)NMFL0_S(树)丰度图 (i)真实(水)丰度图 (j)NMF(水)丰度图 (k)L1/2NMF(水)丰度图 (l)NMFL0_S(水)丰度图图2 基于几种算法下的丰度图像

为了更清晰的比较几种算法解混的精度，本文还将几种算法的端元光谱图像和丰度图像进行了比较分析，如图3～4所示。图3中给出了几种算法(NMF，L1/2NMF，NMFL0_A)对Jasper高光谱数据进行的混合像元分解后的端元估计值与真实端元值的比较图。从图3中可以看出，三种算法的端元光谱基本都接近真实的光谱，尤其在Soil中拟合效果最好。在Tree和Water光谱拟合中，NMF和L1/2NMF算法与真实光谱偏差相对NMFL0_A较大，所以此时NMFL0_A的鲁棒性最好，优势明显。

图3 基于几种算法下的端元信息图

图4是几种算法(NMF，L1/2NMF，NMFL0_S)对Jasper高光谱数据进行的混合像元分解后的端元对应的丰度图像。从图4中可以看出在水部分，三种算法都能较好的接近真实的丰度图。在水和路部分，三种算法差不多，在树部分，NMFL0_S算法的丰度分布图(m)明显比L1/2NMF的丰度图(i)更加接近真实丰度图，在沙土部分相比于NMF算法的丰度分布图(g)，NMFL0_S算法的丰度分布图(o)要清晰很多，更接近真实丰度图。所以综合图3～4，可以看出NMFL0_S和NMFL0_A这两种算法相比NMF和L1/2NMF算法更加优越，效果更好。

(a)真实(树)丰度图 (b)真实(水)丰度图 (c)真实(沙土)丰度图 (d)真实(路)丰度图 (e)NMF(树)丰度图 (f)NMF(水)丰度图 (g)NMF(沙土)丰度图 (h)NMF(路)丰度图 (i)L1/2NMF(树)丰度图 (j)L1/2NMF(水)丰度图 (k)L1/2NMF(沙土)丰度图 (l)L1/2NMF(路)丰度图 (m)NMFL0_S(树)丰度图 (n)NMFL0_S(水)丰度图 (o)NMFL0_S(沙土)丰度图 (p)NMFL0_S(路)丰度图图4 基于几种算法下的丰度图像

4 结束语

本文就NMF的L0约束提出了一种近似NMF的框架(NMFL0_S，NMFL0_A)，它分别约束基础矩阵和系数矩阵，通过推导该算法的稀疏编码方式和矩阵更新的迭代关系，将其与不受约束的NMF技术结合。在高光谱混合像元分解中，实现了对端元矩阵和丰度矩阵的联合稀疏约束，使混合像元分解的精度和效率明显提升，很好的保证了NMF倾向于返回其输入数据的稀疏和基于部分表示的特点。最后采用两组真实仿真数据验证了该算法的优越性和有效性，相对其他算法在光谱解混中的精度和运算速度都具有明显的提升。以上算法是针对L0约束的NMF算法提出的改进算法，如何添加约束项，使稀疏效果更好，运算速度更快的NMF算法是今后有待研究的工作。