基于稀疏和正交约束非负矩阵分解的高光谱解混

陈善学 储成泉

摘 要：针对基于非负矩阵分解（NMF）的高光谱解混存在的容易陷入局部极小值和受初始值影响较大的问题，提出一种稀疏和正交约束相结合的NMF的线性解混算法SONMF。首先，从传统的基于NMF的高光谱线性解混方法出发，分析高光谱数据本身的理化特性;然后，结合丰度的稀疏性和端元的独立性两个方面，将稀疏非负矩阵分解（SNMF）和正交非负矩阵分解（ONMF）两种方法结合应用到高光谱解混当中。模拟数据和真实数据实验表明，相比顶点成分分析法（VCA）、SNMF和ONMF这三种参考解混算法，所提算法提高了线性解混的性能;其中，评价指标光谱角距离（SAD）降低了0.012～0.145。SONMF能够结合两种约束条件的优势，弥补传统基于NMF线性解混方法对高光谱数据表达的不足，取得较好的效果。

关键词：非负矩阵分解; 高光谱解混; 稀疏; 正交;独立性

中图分类号： TN79

文献标志码：A

Hyperspectral unmixing based on sparse and orthogonal constrained non-negative matrix factorization

CHEN Shanxue， CHU Chengquan*

School of Communication and Information Engineering， Chongqing University of Posts and Telecommunications， Chongqing 400065， China

Abstract： Aiming at the problem that hyperspectral unmixing based on Non-negative Matrix Factorization （NMF） is easy to fall into local minimum and greatly affected by initial value， a linear unmixing algorithm based on Sparse and Orthogonal constrained Non-negative Matrix Factorization （SONMF） was proposed. Firstly， based on the traditional NMF hyperspectral linear unmixing method， the physical and chemical properties of the hyperspectral data was analyzed. Then the sparsity of the abundance and the independence of the endmember were combined together， two methods of Sparse Non-negative Matrix Factorization （SNMF） and Orthogonal Non-negative Matrix Factorization （ONMF） were combined and applied into hyperspectral unmixing. The experiments on simulation data and real data show that， compared with the three reference unmixing algorithms of Vertex Component Analysis （VCA）， SNMF and ONMF， the proposed algorithm has improved the performance of linear unmixing， in which the Spectral Angle Distance （SAD） is reduced by 0.012 to 0.145. SONMF can combine the advantages of the two constraints to make up for the lack the expression of hyperspectral data by traditional NMF based linear unmixing methods， and achieve good results.

Key words： Non-negative Matrix Factorization （NMF）; hyperspectral unmixing; sparsity; orthogonality; independence

0 引言

近年來，传感器技术的快速发展大大加快了高光谱遥感图像处理技术在遥感领域的实际应用，高光谱遥感技术的快速发展使得遥感领域的许多方向发生了翻天覆地的变化。特别是在对生活环境特别重视的现在，高光谱技术在探测地表环境、地质灾害等领域发挥着越来越关键的作用[1]。但是，由于遥感传感器的空间分辨率较低，相邻物质可能混合在一起形成混合像元。高光谱解混[2]就是将混合像元分解为各种地物的端元光谱以及对应的丰度系数。由于在混合像元分解时有关端元和丰度的任何信息都是未知的，而且受到复杂的自然条件和地物环境等因素的影响，混合像元分解问题变得更加困难，所以在解决混合像元分解问题之前，首先要确定光谱混合模型。

线性光谱混合模型（Linear Spectral Mixed Model， LSMM）[3]假定端元之间互相独立，是一个描述混合像元现象的有效模型。非负矩阵分解（Non-negative Matrix Factorization， NMF）[4]在基于LSMM的高光谱解混中得到广泛的关注，它采用乘子算法将数据分解为两个非负矩阵，不需要每个端元都存在纯像素，非负的约束自然保证了LSMM中丰度的非负约束（Abundance Non-negativity Constraint， ANC）;然而，由于函数本身的非凸性，NMF具有大量的局部极小值，这将影响最终的解混结果。为了解决这个问题，需要在NMF模型中引入额外的结构性约束。例如，利用高光谱的空间特性[5-6]，将空间信息的约束引入到基于NMF的方法中，或者充分利用高光谱固有的稀疏性特性[7-8]，将稀疏性约束引入到基于NMF的解混方法中。近年来，一种新的稀疏约束NMF （L1/2-NMF）解混方法[9]在高光谱解混中得到了广泛的应用。L1/2-NMF已经被证明是相对较优的方法。然而，与其他稀疏约束NMF一样，L1/2-NMF的解混性能仍然受到初始设置和噪声干扰的影响[10]。

考虑到稀疏约束NMF的高光谱解混方法存在的问题，本文采用经典的顶点成分分析法（Vertex Component Analysis，VCA）[11]得到端元，通过全约束最小二乘法（Fully Constrained Least Square， FCLS）[12]得到丰度，将该算法得到的端元和丰度矩阵作为算法输入的初始值，同时在NMF中加入了正交性约束即正交非负矩阵分解（Orthogonal Non-negative Matrix Factorization， ONMF）[13-14]。通过进行模拟数据和真实数据实验分析，稀疏和正交约束非负矩阵分解（Sparse and Orthogonal constrained Non-negative Matrix Factorization， SONMF）在对高光谱数据进行解混时能够充分利用数据本身的特殊性，在保证端元相互独立同时减小噪声对解混结果的影响。本文的技术路线如图1所示。

1 线性光谱混合模型

在线性光谱混合模型（LSMM）中，每一个端元的光谱特性是其端元光谱特性和丰度比例的线性组合，用矩阵表示为：

R = EA + N  （1）

s.t.  A ≥0， 1 Tp A = 1 Tn

式中： R =[ r i，…， r n]∈ R d×n表示高光谱图像的n个像元以及d个波段; E =[ e 1， e 2，…， e p]∈ R d×p表示混合矩阵包含p个端元; A =[ a 1， a 2，…， a n]∈ R p×n表示端元的丰度矩阵; N 表示误差矩阵。丰度的物理意义明确了丰度矩阵中的每个元素非负，且每一列的和为1。这两个约束条件可以用以下的公式表示：

a i≥0; i=1，2，…，n

（2）

∑ n i=1  a i=1

（3）

2 非负矩阵分解与线性解混的结合

非负矩阵分解的数学模型是将高维的原矩阵分解为两个低维非负矩阵，分解的过程实际上是一个不断迭代优化的过程，建立一个合理的目标函数，再交替迭代求得分解后的非负矩阵，形如：

V ≈ WH

（4）

矩阵 V ∈ R m×n分解成两个矩阵 W 和 H ，这里的 W ∈ R m×t， H ∈ R t×n，并且满足t

非负矩阵分解的表达形式十分契合LSMM，所以可以将非负矩阵分解的方法应用于高光谱解混的条件，由于非负矩阵分解是个NP-hard问题，解决这类优化问题需要引入目标函数，文献[4]中提出了一种基于欧氏距离的目标函数

J（ W ， H ）= 1 2 ‖ Y - WH ‖2F

（5）

其中，‖·‖2F表示Frobenius范数，因为目标函数是非凸的，所以寻找全局的最优解很困难的。因此，使用迭代方法如乘法更新规则被认为是解决这类基于NMF的问题的有效工具之一。式（5）的乘法更新法则可以表示为：

W ← W .*  VH T  WHH T  H ← H .*  W T V   W T WH

（6）

2.1 基于稀疏约束的非负矩阵分解

在高光谱解混中，其中经常被用到的约束条件是基于丰度稀疏性的约束，稀疏性约束指的是每个混合像元中包含的端元数量远小于总的端元数量，在传统的稀疏约束的NMF（Sparse Non-negative Matrix Factorization， SNMF）解混方法中，L1 / 2正则化[9]作为稀疏约束项被引入NMF以追求更稀疏的丰度表示，由于L1/2-NMF能够有效利用数据本身的稀疏性，与其他稀疏NMF方法相比，它显示出较大的优势，基于L1/2-NMF的目标函数可以表示为：

J（ W ， H ）= 1 2 ‖ Y - WH ‖2F+α‖ H ‖1/2

（7）

式中：‖ H ‖1/2=∑ L，N l，n=1 （ h l×n）1/2 表示正则化，α∈ R +表示正则化参数。

L1/2-NMF算法解混的效果受到端元以及丰度初始设置值的影响较大，另外，噪声对解混结果也有很大的影响，所以为了取得更好解混效果，需要将更优的约束条件引入到NMF目标函数中去。

2.2 基于稀疏和正交约束的非负矩阵分解

由于高光谱解混提取的端元是相互独立的，鉴于这样的特殊性，本文将正交约束[13]应用到非负矩阵分解中，该目标函数可以表示为：

J（ W ， H ）= 1 2 ‖ Y - WH ‖2F （8）

s.t.  W ≥0， H ≥0， WW T= I

正交约束被证明等价于K-mean聚类等价，式（8）中对端元矩阵 W 加入正交约束，对提取出来的端元而言，保证了独立性。考虑到高光谱数据集高维度的特点，通常基于正交非负矩阵分解算法复杂度为O（m2n），每当m增加时，相应的算法也会变得更加复杂，直接引入正交约束之后需要一定的先验知识并且需要很大计算量，迭代求解的复杂度也随之增加，所以将正交非负矩阵分解应用到高光谱解混中时需要降低算法复杂度。将正交约束直接表示在目标函数中，将正交约束作为目标函数的一部分，不需要借助额外的约束条件，这种方法可以一定程度上减少计算量，改进后的基于正交约束的NMF的目标函数可以表示为：

J（ W ， H ）= 1 2 ‖ Y - WH ‖2F+β‖ W T W - I ‖22 （9）

s.t.  W ≥0， H ≥0

式中：β表示正交回歸参数，‖ W T W - I ‖22表示2范数约束。

通過对式（5）引入端元正交性约束，可将丰度稀疏性和端元正交性约束非负矩阵分解（SONMF）的目标函数表示为：

J（ W ， H ）=  1 2 ‖ Y - WH ‖2F+α‖ H ‖1/2+

β‖ W T W - I ‖22

（10）

s.t.  W ≥0， H ≥0， 1 Tp H = 1 Tn

对于式（10）的求解，这里使用乘法更新法则对矩阵 W 和 H 进行更新。首先对目标函数分别求关于 W 和 H 的偏导数，因为矩阵 W 和 H 是非负的，记

ξ W =  W   WHH T+β WW T W

ξ H =  H   W T WH + α 2  H - 1 2

（11）

将ξ W ，ξ H 代入：

W = W -ξ W J W

H = H -ξ H J H

（12）

使用最小二乘迭代规则，得到 W 和 H 的更新规则：

W ← W .*  VH T+2β W   WHH T+2β WW T W

H ← H .*  W T V   W T WH + α 2  H - 1 2

（13）

2.3 预处理环节

将SONMF算法应用于高光谱解混时，需要考虑一些预处理的环节。首先是端元数目的确定，本文采用最小误差信号子空间识别（Hyperspectral Signal Identification by Minimum Error，HySime）[15]方法确定端元的数目;接着是端元和丰度初始化的问题，考虑到基于范数约束NMF受初始值影响较大，

用顶点成分分析法得到端元矩阵，通过全约束最小二乘法得到丰度矩阵，将得到的端元和丰度矩阵作为算法输入的初始值，这样不仅可以降低初始值对解混的影响，还可以加快算法的执行速度，减少算法的执行时间。关于参数设置的问题，稀疏正则化参数α主要根据矩阵丰度稀疏度来设置，具体的表达式如下：

λ= 1  L  ∑ l   N -‖ x l‖1/‖ x l‖2  N -1

（14）

其中 x l表示高光谱数据第l维的光谱。正交化参数的设置根据文献[14]设为0.05。

最后是停止条件设置的问题，本文停止算法的条件是迭代次数达到500时停止。

SONMF算法概括如下：

1）输入图像矩阵 Y ，用HySime算法估计端元的个数，采用VCA初始化端元矩阵 W ，利用FCLS初始化丰度矩阵 H ，设定参数α，β。

2）将初始化的矩阵 W 和 H 代入乘法更新式（13）进行更新。

3）如果目标函数式达到最大的迭代次数Tmax，算法停止运行，否则返回2）。

2.4 收敛性的证明

正交约束目标函数式（9）收敛性的证明参照了文献[4]的思路，证明目标函数不递增，考虑矩阵 W 中的某个元素wij，对式（9）求 J 对wij的一阶偏导数和二阶偏导数可得：

J′ij（w）= WHH T+2β WW T W - YH T-2β W

（15）

J″ij（w）=（ HH T）ij+2β（（ W T W ）jj+w2ii-1）

（16）

构造关于 J 的辅助函数，记

F（w，wtij）=Jij（wtij）+J′ij（wtij）（w-wtij）+

（ WHH T）ij+2β（ WW T W ）ij wtij  （w-wtij）2

（17）

证明收敛性只要证明F（w，wtij）≥Jij（w）即可。

F（w，wtij）=Jij（wtij）+J′ij（wtij）（w-wtij）+

（ WHH T）ij+2β（ WW T W ）ij wtij  （w-wtij）2

（18）

根据条件（ W T W ）jj=1， Jij（w） 可以表示为：

Jij（w）= Jij（wtij）+J′ij（w-wtij）+

[（ HH T）ij+2βw2ii]（w-wtij）2

（19）

同时可以得到：

（ WHH T）ij=∑ k p=1 wtip（ HH T）pj≥wtij（ HH T）jj

（20）

（ WW TW）ij=∑ k p=1 wip（ W T W ）pj≥wii（ W T W ）≥     wii （ ∑ k p=1 wpiwpj ） ≥wiiwiiwij=w2iiwij

（21）

通过比较式（17）和（19）得到F（w，wtij）≥Jij（w），所以原目标函数式是收敛的。

3 实验分析

3.1 精度评价指标

衡量高光谱解混的精度通常通过两种评价指标：光谱角距离（Spectral Angle Distance，SAD）和均方根误差（Root Mean Square Error，RMSE）。

SAD的计算公式如下：

RSAD（ a ， b ）=arccos   a T b  ‖ a ‖2·‖ b ‖2

（22）

其中： a 表示提取的端元光譜矢量， b 表示参考光谱矢量即光谱库中的已知光谱矢量。

RMSE的计算公式如下：

RRMSE= 1 n ∑ n i=1   1 L ‖ r i-   i‖22

（23）

其中：{   i}ni=1表示重构的图像。RMSE的值越小表示原始图像与重构图像间的误差越小，解混的精度越高。

3.2 合成数据分析

本实验在处理器i5，CPU频率3.5GHz，内存8GB的计算机上使用Matlab2014a操作平台进行实验，合成数据是从美国地质勘探局USGS光谱中选取五个端元Alunite、Buddingtonit、Calcite、Kaolinite和Muscovite，它们的光谱包含224个波段，波长范围为0.38～2.5μm，光谱图见图2。

将这五个端元以Dirichlet分布的形式进行混合，同时对端元丰度的和进行归一化操作，为了模拟数据采集过程，按照指定的信噪比（Signal-to-Noise Ratio，SNR）在模拟图像中加入不同的噪声等级，其中SNR可表示为：

SNR=10 lg E[ x T x ] E[ n T n ]

（24）

式中： x 和 n 表示原始信号和相应的噪声，E[·]表示期望值。

实验比较了VCA[11]、ONMF[13]、SNMF[9]以及本文采用的SONMF算法的性能。

实验1

不同噪声条件下算法鲁棒性分析。本实验研究高斯噪声对实验的影响，所以保持图像大小为100×100，端元数目为5，SNR从15 dB到35dB，间隔为5dB。

实验2

不同像元数目条件下对实验的影响。本实验研究像元数目对实验的影响，所以保持端元数目为5，SNR为20dB的条件不变，图像大小分别为40×40，60×60，80×80，100×100。

实验结果如图3、4所示，可以看出，SONMF解混效果整体要优于SNMF和ONMF的解混效果。

本文采用的数据集是1997年机载可见光及红外光谱成像仪采集到的Urban地区高光谱图像（307×307像素），总共有210个波段。除去低信噪比和水汽吸收的波段（1-4，76，87，101-111，136-153，198-210），只有162个有效波段，该地区包含了屋顶、树木、草地、道路等多种地物的混合，图5显示了四种地物光谱的真实丰度图。

首先通过HySime估计端元数目为4，再使用VCA算法提取的端元及FCLS估计的丰度作为初始值，分别使用SNMF、ONMF和SONMF算法对该地区的高光谱图像进行处理，表1给出了不同算法最终的SAD的标准值和均值，加粗的字体是提取较好的结果，本文提出的SONMF算法的解混结果是比较理想的，均值SAD降低0.012～0.145。

图6给出了这四种算法对Urban地区高光谱数据进行解混得到的端元对应的丰度图。

将图6中提取的光谱丰度图与图5中地物光谱的真实端元丰度图进行对照可以看出，除了VCA算法分解出的端元丰度图效果较差之外，其他三种算法分解出的端元丰度图分布在道路和草地部分与真实端元丰度图分布比较接近;然而在屋顶和树木部分，SONMF算法分解出端元丰度图明显要优于其他的算法。

4 结语

混合像元处理是高光谱遥感领域的一个热门研究，本文采用线性光谱混合模型来解决光谱分解的问题，重点讨论了非负矩阵分解的方法，非负矩阵分解没有独立的解决方案，需要在目标函数中结合其他的约束条件求解，本文将稀疏性与端元正交性约束结合作为非负矩阵分解的约束条件，同时在迭代初始值方面作了一些改进。本文算法在保证提取端元的独立性和稀疏端元本身丰度的同时，减小了噪声和局部极小值对解混解混结果的影响，取得了良好的解混效果。但本文算法相对比较复杂，如何进一步通过预处理的手段降低算法复杂度是下一步的研究工作。

参考文献

[1] 童庆禧，张兵，郑兰芬.高光谱遥感——原理、技术与应用[M].北京：高等教育出版社， 2006：15-30. （TONG Q X， ZAHNG B， ZHENG L F. Hyperspectral Remote Sensing——Principle， Technology and Application [M]. Beijing： Higher Education Press， 2006： 15-30.）

[2] BIOUCAS-DIAS J M， PLAZA A， DOBIGEON N， et al. Hyperspectral unmixing overview： geometrical， statistical， and sparse regression-based approaches [C]// Proceedings of the 2011 IEEE International Geoscience and Remote Sensing Symposium. Piscataway， NJ： IEEE， 2011： 1135-1138.

[3] 张兵，孙旭.高光谱图像混合像元分解[M].北京：科学出版社， 2015：6-22. （ZAHNG B， SUN X. Hyperspectral Image Unmixing [M]. Beijing： Science Press， 2015： 6-22.）

[4]  LEE D D， SEUNG S S. Learning the parts of objects by non-negative matrix factorization [J]. Nature， 1999， 401： 788-791.

[5] 袁博.馬尔可夫随机场的空间相关模型在非负矩阵分解线性解混中的应用[J].计算机应用，2017，37（12）：3563-3568. （YUAN B. Application of MRFs spatial correlation model in NMF-based linear unmixing [J]. Journal of Computer Applications， 2017， 37（12）： 3563-3568.）

[6] WANG X， ZHONG Y， ZHANG L， et al. Spatial group sparsity regularized nonnegative matrix factorization for hyperspectral unmixing [J]. IEEE Transactions on Geoscience and Remote Sensing， 2017， 55（11）：6287-6304.

[7] RAJABI R， GHASSEMIAN H. Spectral unmixing of hyperspectral imagery using multilayer NMF [J]. IEEE Geoscience and Remote Sensing Letters， 2014， 12（1）： 38-42.

[8] ZHAO Y， ZHOU Z， WANG D， et al. Hyperspectral image unmixing algorithm based on endmember-constrained nonnegative matrix factorization [J]. Frontiers of Optoelectronics， 2016， 9（4）： 627-632.

[9]  QIAN Y， JIA S， ZHOU J， et al. Hyperspectral unmixing via L1/2 sparsity-constrained nonnegative matrix factorization [J]. IEEE Transactions on Geoscience & Remote Sensing， 2011， 49（11）： 4282-4297.

[10] 余肖玲.非负矩阵分解理论及其在高光谱解混中的应用[D]. 成都：成都理工大学， 2015：45-60. （YU X L. The theory of Nonnegative matrix factorization and its application in the hyperspectral unmixing [D]. Chengdu： Chengdu University of Technology， 2015： 45-60.）

[11] NASCIMENTO J M P，BIOOUCASDIAS J M B. Vertex component analysis： a fast algorithm to unmix hyperspectral data [J]. IEEE Transactions on Geoscience and Remote Sensing， 2005， 43（4）： 898-910.

[12] HEINZ D C， CHANG C. Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery [J]. IEEE Transactions on Geoscience and Remote Sensing， 2002， 39（3）： 529-545.

[13] POMPILI F， GILLIS N， ABSIL P， et al. Two algorithms for orthogonal nonnegative matrix factorization with application to clustering [J]. Neurocomputing， 2014， 141：15-25.

[14] 李孟杰，谢强，丁秋林.基于正交非负矩阵分解的K-means聚类算法研究[J].计算机科学，2016，43（5）：204-208. （LI M J， XIE Q， DING Q L. Orthogonal Non-negative Matrix Factorization for K-means Clustering [J]. Computer Science， 2016， 43（5）： 204-208.）

[15] BIOUCASDIAS J M， NASCIMENTO J M P. Hyperspectral subspace identification [J]. IEEE Transactions on Geoscience and Remote Sensing， 2008， 46（8）：2435-2445.