非负矩阵分解与光谱解混

孙莉,于瑞林,吴杰芳

非负矩阵分解与光谱解混

孙莉1,于瑞林1,吴杰芳2

1. 山东农业大学 信息科学与工程学院, 山东 泰安 271018 2. 泰山学院 数学与统计学院, 山东 泰安 271000

非负矩阵分解(NMF)用两个非负矩阵的乘积近似原始数据对应的非负矩阵，它为基于线性光谱混合模型的光谱解混提供了新途径。给出NMF在光谱解混中三个矩阵的具体含义后，用五种求解NMF的有效算法，对Jasper Ridge的高光谱遥感图像进行解混。讨论了五种算法的迭代方式以及收敛性质。实验结果表明，五种算法能成功分离出4种端元光谱以及相应的丰度谱图，其中有效集型算法表现突出。

非负矩阵分解; 光谱解混; 界约束优化; 有效集

非负矩阵分解(Nonnegative matrix factorization, NMF)的思想可以追溯到1994年Pattero和Tapper的文章[1]，随后1999年，Lee DD和Seung HS独立的提出了NMF的概念后[2]，该方法迅速引起了众多学者的关注，蓬勃发展起来。用两个非负矩阵的乘积近似原始数据对应的非负矩阵，非负性的要求保证了通过部分表示原始数据的形式，即样本数据只允许加性组合，这使得分解结果对原始数据的特征及结构具有相当的表达能力，有利于发现数据的局部特征。

1 NMF的基本思想与主要求解算法

1.1 NMF的基本思想

若的各列表示每一样本的相关属性，则经NMF分解得到的矩阵称为基矩阵，称为系数矩阵。将和按列分块，即=(1,2,…,w)，=(1,2,…,v)，由(1)可得个等式，

w≈Uv，=1,2,…,(2)

(2)式说明，每一个样本w可近似的看作非负矩阵的列向量的非负线性组合，组合系数是v的分量。于是，矩阵可以看作是对原始矩阵进行线性逼近的一组基，而则是样本集在基上的非负线性投影系数。寻找合适的基向量组，用其表征大量的数据向量，一方面可以给出数据之间潜在的结构关系，另一方面实现了数据压缩，用少量的数据实现对原始数据很好的逼近。通过NMF对矩阵进行分解后，所得矩阵和在不同的应用领域有不同的含义，它在图像分析、文本聚类、语音处理、生物医学等领域应用广泛。本文第二部分，我们会给出NMF应用于光谱解混时，和的具体含义。

1.2 NMF问题模型

求解矩阵和的常用方法是通过最小化和的差距，即

1.3 求解NMF问题的有效算法

针对NMF问题模型，已经提出了多种行之有效的求解方法，现有算法可分为两类：第一类，基于梯度下降的交替迭代公式，通过公式交替实现和的更新。目前求解NMF的常用方法——乘法迭代公式便属于这一类。算法的特点是结构简单，易于实现，但其在运行效率与收敛性分析上有所欠缺。第二类，交替非负最小二乘(Alternating Nonnegative Least Squares, ANLS)框架下的有效算法，ANLS算法框架如下：

ANLS算法框架

步骤 1 给定初始矩阵0和0，:=0。

步骤 2 依次实现和的更新如下:

基于ANLS算法框架下的优化方法有投影梯度法、投影(拟)牛顿法、有效集型方法和块转轴方法等。2007年Lin CJ的文献[3]是利用传统界约束优化算法求解NMF问题的开端。以优化理论为基础，这些算法具有良好的收敛性质。

下面给出5种求解NMF的算法。

（1）乘法迭代公式(MU)

该方法一经提出，便成为求解NMF问题的标准算法[4]，其迭代规则为：

（2）交替最小二乘算法(ALS)

通过交替求解最小二乘问题实现和的更新，一旦矩阵元素不满足非负性要求，直接将该元素定义为0。忽略非负性约束，无法保证目标函数的单调下降，因此无法获得算法的收敛性结论，但这一算法运行速度快，其迭代规则如下：

和已知，通过UUV=UW，求出，将中的负数直接赋值为0；

和已知，通过VVU=VW，求出，将中的负数直接赋值为0。

（3）投影梯度法 (PGM)

（4）投影牛顿法 (PNM)

Gong PH和Zhang CS以Lin CJ的投影梯度法为基础给出了投影牛顿法[5]，它通过乘子信息同样将变量分为有效变量和非有效变量，以投影牛顿方向为搜索方向实现非有效变量的更新。这一方法的优点在于一步迭代允许有多个指标进出有效集，且搜索方向的定义借助于目标函数的二阶信息，它获得了二阶收敛的速度。

（5）有效集型算法 (AS)

Kim H和Park H提出了求解NMF的有效集型算法[6]，他们将待求变量分为两类：有效变量(非负约束中等式成立的变量)和非有效变量，对非有效变量通过求解无约束的最小二乘问题实现更新。随后以目标函数下降为前提，重新制定有效变量和非有效变量。算法是否停止迭代，通过检测最优化条件判断，即非有效变量的非负性和有效变量对应乘子的非负性是否满足。这一算法的主要问题在于每步迭代仅允许一个指标进出有效集，若初始矩阵的选取不合适，算法的计算复杂度很高。

2 光谱解混中的NMF

2.1 遥感图像的矩阵表示

高光谱图像可以看成是三维数据构成的立方体，它的和维表现的是地球表面的坐标，第三维是波段，它由光波的频率所决定。例如，将一幅AVIRIS图像读取为××的矩阵，其中×对应图像的宽度和高度，对应波段数，于是一幅AVIRIS图像对应了幅×的二维图像，每幅图像对应地球上相同位置在不同波长下的影像。不妨设图像的像素点有个，=×，于是将上述矩阵重新排列可得到×的矩阵。的行向量代表某一像素点在不同波段的灰度值(DN值)，而的列向量对应某一波段下，个像素点的DN值。

2.2 遥感图像的矩阵表示

通常遥感影像以像元为单位获取地物信息，所记录的是像元内所有地物光谱的综合。在实际中，由于空间分辨率的限制以及地物的复杂多样性，一个像元内往往会包含多种地物类型，称为混合像元。为提高分类结果对真实地表覆盖的描述准确性，需要进入像元内部，将混合像元分解为不同的“基本组成单元”或“端元”，并求得这些基本组分所占的比例，这就是所谓的“光谱解混”过程。

在线性光谱混合模型中，混合光谱是端元光谱与其比例的线性组合，其数学表达式如下：

其中，表示端元的个数，为维光谱向量对应图像中某一像元在个波段的DN值，1,2,…,e分别对应个端元对应的光谱向量，c代表像元中端元e所占的比例，为误差项。

2.3 遥感图像的矩阵表示

非负矩阵分解与光谱解混中线性光谱混合模型相吻合，为光谱解混提供了一种新的途径。利用NMF实现的分解后，得到的矩阵为丰度矩阵，它的行给出某一像元中各端元所占比例，列向量则对应某个端元的空间分布。为端元矩阵，它的行向量为某个端元的光谱向量[7,8]。

3 数值比较

本文所采用的数据来自于2006年机载可见光及红外成像光谱仪(AVIRIS)采集到的Jasper Ridge高光谱图像，图像大小为512×614像素，自380 nm至2500 nm，有224个波段，光谱分辨率达到9.46 nm。我们截取了100×100像素的子图，该区域包含了树木、道路、裸土等多种地物的混合，去除信噪比低和水蒸气吸收波段(1~3，108~112，154~166以及220~224)，余下198个有效波段。图1为第35，20，7波段图像合成的伪彩色图像。

图1 Jasper Ridge的AVIRIS伪彩色图

Jasper Ridge地区的树木、水体、裸土以及道路，4种地物的真实丰度谱图，如图2所示。

图2 Jasper Ridge真实地物丰度谱图

设置端元数目为4，即=4，利用NMF实现上述遥感图像的光谱解混，算法程序由Matlab 7.0编写。我们主要从解混后的丰度谱图，分解所得的端元光谱与真实端元光谱的余弦值以及运行时间3个方面比较5种算法的性能。求解NMF的5种算法分别为乘法迭代公式(MU)，交替最小二乘算法(ALS)，投影梯度法(PGM)，投影牛顿法(PNM)以及有效集型算法(AS)。

NMF问题(2)非凸，导致算法的运行效率与初始点的选取密切相关。数值测试中，利用(,)，(,)随机生成初始矩阵0和0，随后5种算法以相同的0，0开始运行。算法的最高迭代次数设置为1000，运行时间上限为104s，PGM和PNM增加了下述终止准则，即：

下标表示的第个分量。

图3 采用5种算法获得的Jasper Ridge四地物丰度谱图比较

图3为采用5种算法对图像数据矩阵进行NMF分解，所得丰度矩阵对应的丰度谱图，图中第1至4列分别对应树木、水体、裸土以及道路。通过对比发现，5种算法均成功获得树木与裸土的丰度谱图。其中MU，PGM，PNM，AS所得树木的丰度谱图与真实地物吻合度较高，同时，裸土丰度谱图以ALS和AS算法所得为佳。水体与道路的分离不显著，如图3所示，第二列图像对应的水体丰度谱图中，水体与道路的信息一并显示出来。通过NMF分解所得的端元矩阵，提供了四种端元：树木，水体，裸土以及道路的端元光谱信息。下面计算分离的光谱信息与美国地质调查局(USGS)库中的参考光谱间的光谱角距离。光谱角距离(Spectral Angel Distance, SAD)，即两个光谱向量之间的夹角，计算公式如下：

其中V，V表示不同端元的光谱向量，通过SAD计算出的角度越小代表与参考光谱越匹配。

表4 USGS谱库的真实光谱和提取的端元光谱的SAD比较

由表1看出AS算法提取出的4个端元光谱与真实光谱最接近，PGM和PNM算法所得的丰度谱图以及4个端元光谱信息十分接近。事实上，PGM和PNM算法的终止准则以及算法框架大致相同，主要区别在于搜索方向的定义，PGM借助于投影梯度方向，而PNM则建立在投影牛顿方向的基础上。最后，利用5种算法完成光谱解混工作所需的计算时间如表2所示，时间单位：秒。

表5 不同算法完成光谱解混的时间比较

PNM采用了借助二阶信息的搜索方向，计算时间较之PGM明显减少。5种算法中AS算法的运行时间最少，经分析，AS算法在两种情形会终止迭代，一是超出最大迭代次数；二是满足一阶最优性条件，即非有效变量非负性以及有效变量对应乘子的非负性满足。本测试在第二种情形终止迭代。

4 小结

NMF应用于光谱解混工作，在给定端元个数的条件下，可同时提取出端元光谱和每个端元的丰度谱图。本文总结了5种求解NMF的有效方法，并具体应用于高光谱图像的解混。结果表明，有效集型算法的分解结果与真实地物的吻合度较高，且运行时间有明显优势。我们将进一步将有效集识别方法[9]推广至正则化NMF问题[7,8]的求解。

[1] Paatero P, Tapper U. Positive matrix factorization: a non-negative factor model with optimal utilization of error estimates of data values[J]. Environmetricsa, 1994,5(2):111-126

[2] Lee DD, Seung HS. Learning the parts of objects by non-negative matrix factorization[J]. Nature, 1999,401(6755):788-791

[3] Lin CJ. Projected gradient methods for non negative matrix factorization[J]. Neural Comput.,2007,19(10):2756-2779

[4] Lee DD, Seung HS. Algorithms for non-negative matrix factorization[C]//Neural Information Processing Systems Conference. Massachusetts USA: MIT Press, 2001:556-562

[5] Gong PH, Zhang CS. Efficient nonnegative matrix factorization via projected Newton method[J]. Pattern Recognition, 2012,45(9):3557-3565

[6] Kim H, Park H. Nonnegative matrix factorization based on alternating nonnegativity constrained least squares and active set method[J]. SIAM Journal on Matrix Analysis and Applications, 2008,30(2):713-730

[7] Gillis N, Plemmons RJ. Sparse nonnegative matrix underapproximation and its application to hyperspectral image analysis[J]. Linear Algebra and its Applications, 2013,438(10):3991-4007

[8] 孙莉,赵庚星.基于正交非负矩阵分解的高光谱遥感图像混合像元分解[J].山东省农业大学学报(自然科学版),2017,48(2):264-267

[9] 孙莉,贺国平,房亮.基于求解大规模界约束问题的三种有效集识别策略的比较[J].数值计算与计算机应用,2009,30(1):41-47

Non-negative Matrix Factorization and Hyperspectral Unmixing

SUN Li1, YU Rui-lin1, WU Jie-fang2

1.271018,2.271000,

Nonnegative matrix factorization approximates the data matrix with the production of two low-rank matrices. It provides a new approach to hyperspectral image analysis, which is based on the standard linear mixing model. After introducing the implications of the three matrices in NMF, we summarize five different computational algorithms which solve NMF successfully, and their iteration procedures and convergence properties are analyzed. In the experimental test, five algorithms are employed to unmix the hyperspectral image Jasper Ridge. Numerical results show that, they can detect the endmembers spectral signatures and the corresponding fractional abundance. A good choice is the active set type method which performs noticeably well.

Nonnegative matrix factorization; hyperspectral unmixing; bound constrained optimization; active set

O221;TP751.1

A

1000-2324(2019)05-0908-05

10.3969/j.issn.1000-2324.2019.05.038

2015-09-05

2015-10-29

国家自然科学基金项目(11701337,10901094,11301307);山东省自然科学基金项目(ZR2016DM03);山东省高等学校科技计划项目(J16LI16);大学生创新创业训练计划项目(201410434097)

孙莉(1980-),女,博士,副教授,研究方向为最优化算法与理论. E-mail:sunlishi@hotmail.com