光子学报  2018, Vol. 47 Issue (3): 0310002  DOI: 10.3788/gzxb20184703.0310002
0

引用本文  

王瀛, 何欣, 左方. 基于最大整体包容度约束非负矩阵分解的高光谱遥感图像混合像元分析算法[J]. 光子学报, 2018, 47(3): 0310002. DOI: 10.3788/gzxb20184703.0310002.
WANG Ying, HE Xin, ZUO Fang. Mixed Data Analysis Algorithm Based on Maximum Overall Coverage Constraint Nonnegative Matrix Factorization for Hyperspectral Image[J]. Acta Photonica Sinica, 2018, 47(3): 0310002. DOI: 10.3788/gzxb20184703.0310002.

基金项目

国家自然科学基金(No.61703141),河南省基础与前沿技术研究(No.162300410198)和河南省博士后科学基金(No.153040)资助

第一作者

王瀛(1976-), 男, 副教授, 博士, 主要研究方向为高光谱遥感图像处理. Email:wangying@henu.edu.cn

通讯作者

何欣(1974-), 男, 教授, 博士, 主要研究方向为无线传感网与分布式云计算. Email:hexin@henu.edu.cn

文章历史

收稿日期:2017-10-18
录用日期:2017-12-05
基于最大整体包容度约束非负矩阵分解的高光谱遥感图像混合像元分析算法
王瀛1, 何欣2, 左方1    
(1 河南大学 智能网络系统研究所, 河南 开封 475001)
(2 河南大学 河南省现代网络技术实验教学示范中心, 河南 开封 475001)
摘要:针对高光谱遥感图像中存在高度混合无纯像元的现象,提出了端元整体包容度约束,并将其加入非负矩阵分解的目标函数.在满足端元非负性与和为一约束的同时,利用数据在特征空间的几何特性,要求端元构成的单形体所容纳的像元尽可能多.该算法不需对原始数据降维,不损害数据的物理意义,在迭代过程中使用乘性规则,避免了传统梯度优化过程中常见的整体步长难以控制现象.对模拟图像和真实图像进行实验评测并比较了提取端元精准度、鲁棒性以及执行效率,结果表明,本文算法可有效分析高光谱遥感图像混合像元.
关键词高光谱图像    端元    非负矩阵分解    凸面几何学    单形体    
中图分类号:TP751.1      文献标识码:A      文章编号:1004-4213(2018)03-0310002-9
Mixed Data Analysis Algorithm Based on Maximum Overall Coverage Constraint Nonnegative Matrix Factorization for Hyperspectral Image
WANG Ying1, HE Xin2, ZUO Fang1    
(1 Institute of Intelligence Networks System, Henan University, Kaifeng, Henan 475001, China)
(2 Henan Experimental Teaching Demonstration Centre of Modern Network Technology, Henan University, Kaifeng, Henan 475001, China)
Foundation item: The National Natural Science Foundation of China (No. 61703141), the Basic and Frontier Technology Research of Henan Province Science and Technology Department (No. 162300410198) and the Henan Postdoctoral Science Found (No. 153040)
Abstract: In order to analyze hyperspectral images consisted of highly mixed pixels, a new endmembers overall coverage constraint was proposed and introduced in objective function of nonnegative matrix factorization, which forcely maximizes the number of pixels contained in the simplex constructed by endmembers using data geometrical properties in the feature space while satisfies data nonnegative and abundance sum-to-one constraint simultaneously. In the maximum overall coverage constraint nonnegative matrix factorization algorithm, the dimensionality reduction process is prevented to preserve the physical meaning of the source image and multiplicative update rules are applied to avoid stepsize selection problem occurred in traditional gradient-based optimization algorithm frequently. To evaluate the accuracy of endmembers extraction, the performance and robustness, experiments are designed on synthetic and real images. The results demonstrate that the proposed algorithm is an effective method to analyze mixed data in hyperspectral image.
Key words: Hyperspectral image    Endmember    Nonnegative matrix factorization    Convex geometry    Simplex    
OCIS Codes: 100.4145;110.4234;300.6170;100.2960;110.3010
0 引言

高光谱遥感能在空间成像的同时以较高的光谱分辨率记录下成百个连续的光谱通道,将近似连续的地物光谱曲线和地表空间图像结合在一起,所形成的图像立方体信息丰富,可以同时提供反应地物之间空间几何关系的图像信息和用于鉴别单个物质的光谱辐射信息,因而具备广阔的应用前景[1].据其特性,高光谱遥感图像分析工作可在光谱空间、图像空间以及特征空间中展开,特别是在高维特征空间中,一些数学理论和方法得以应用[2].在高光谱图像中,由单一地物构成的像元称为纯像元,亦称端元.研究表明,混合像元和端元之间的关系大致可用线性光谱混合模型描述,寻找端元以及求得混合像元相应的线性表出分别对应着端元提取和丰度反演分析工作[3-4].像元分析的结果是许多后续应用的基础,由于光谱成像仪空间分辨率所限,以及成像区域复杂地物分布的客观现实,像元多以高度混合的情况存在于图像中,复杂混合状态下的像元分析工作是高光谱遥感图像处理中的热门研究方向之一[5].

以往的研究中,以线性混合模型为基础,考虑凸面几何学特性,对光谱数据云团在高维特征空间中的分布加以处理,是做像元分析的一个重要思路.其中涉及到对数据做变换、投影,同时考虑端元处在数据云团边界且混合像元总在端元围成的高维单形体中的事实,提出了一些端元提取算法,包括像元纯度索引(Pixel Purity Index, PPI)[6]及其改进算法快速迭代的像元纯度索引(Fast Iterative Pixel Purity Index, FIPPI)[7]、N-FINDR[8-9]、最小体积变换(Minimum Volume Transform, MVT)[10]、顶点成分分析(Vertex Component Analysis, VCA)[11]等.上述算法在图像中有纯像元存在的前提下,能够快速得到端元集合,然而在高度混合的场合下,这些算法的解并不唯一,也缺乏一个再度求精的依据和步骤,而且不能避免负值的出现.大部分算法因为使用到投影变换手段,对像元的物理意义也有损害.

非负矩阵分解(Nonnegative Matrix Factorization,NMF)在矩阵所有元素均非负的限制下完成对高阶矩阵的低阶逼近.在信息处理领域中很多数据具备非负性的特点,而矩阵的低阶逼近能够降低维数,展现特征[12].高光谱遥感数据同样具备非负性,并且其线性混合模型与非负矩阵分解高度契合,因而近年来发展了一批基于非负矩阵分解的端元提取算法.非负矩阵分解算法是在非负的限制下,通过最小化损失函数完成矩阵分解,达到分解后的整体最优.由于其目标损失函数是非凸的,导致具有不唯一的局部最优解,因此实际在应用非负矩阵分解时,通常需要施加新的约束使得结果更为精准.具体到高光谱图像混合像元分析方面,在保证数据非负性的同时需要考虑和为一约束,以及兼顾端元向量的其他特性.为了应用非负矩阵分解思想完成混合像元分析,Miao和Qi提出了最小体积约束下的NMF(Minimum Volume Constraint Nonnegative Matrix Factorization, MVC-NMF)方法[13],将降维后数据中端元所围成的单形体体积最小作为额外的约束,得到了不错的效果;刘建军等优化了MVC-NMF算法,提出基于交替方向乘子法的非凸项约束NMF迭代算法,加快了求解速度并增加了算法的稳定性[14];耿修瑞等提出的几何优化模型(Geometric Optimization Model,GOM)[15]要求在降维后的低维空间中,任一混合像元和端元围成的单形体体积和最小,其目标函数中仅包含混合像元矩阵;计璐艳引入了像元质心坐标的概念,在做非负矩阵分解的同时,通过计算体积比完成丰度矩阵的估计,提出了BC-NMF算法[16],能够显著地降低计算量;杨斌结合了粒子群算法和最小端元距离约束,提出的约束NMF框架下高维自适应粒子群端元提取算法[17]在一定程度上避免了传统的NMF算法容易陷入局部最优的问题.通过对现有NMF框架下像元分析算法的回顾可以看出,此类算法的研究重点大致集中在两点:一是引入能够更好反映高光谱图像中端元特性的约束项,作为附加惩罚因子加入到目标函数中;二是在传统梯度下降法的思想下,选用合适的迭代步骤保持收敛速度和元素的非负性.现有算法所提出的约束项能够对端元特性不同方面加以反映,但大多数是在降维后的数据上加以处理,这样做的弊端在于降维产生负值的同时改变了像元的物理特性,也改变了像元在特征空间的表现,特别是在像元高度混合的情况下,所输出的端元对原图的表达存在较大误差;另一方面,由于引入额外的约束项,目标函数最小值求解过程中仅依赖整体梯度信息,步长选择缺乏指点,同时非负性维持需要一个强制将负值变为零的过程,最终造成输出结果不确定,乃至和最优解相去甚远.

对于高度混合的局面,端元的整体包容性值得关注,本文在NMF框架下提出了一种代表端元整体包容程度的约束项,反映高度混合图像中端元的整体覆盖性,并且在迭代求解过程中,采用经典的乘法更新规则而非单一的步长控制,以对收敛速度和非负性给予保障.

1 模型和方法 1.1 线性光谱混合模型

线性光谱混合模型指出,高光谱图像中的每个像元都是空间相邻端元的线性组合,并且满足组合丰度系数和为一约束以及所有元素非负约束,即

$ \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{EC}} + \mathit{\boldsymbol{N}} $ (1)

式中,Xn×m维像元矩阵,n为光谱波段数目,m为像元总数,En×p维像元矩阵,其每一个列向量代表一个端元光谱,p为端元个数;Cp×m维丰度矩阵,其列向量为相应像元的端元组合丰度;N为误差.所有元素均非负并对C中每一个列向量cj, 1≤jm有如下约束

$ \sum\nolimits_{i = 1}^p {{\mathit{\boldsymbol{c}}_{i,j}} = 1} $ (2)

满足上述模型的像元数据在特征空间中处于以端元为顶点的单形体中,这是利用凸面几何学思想求解端元问题的理论基础.像元分析需要完成矩阵EC的求解,在纯像元不存在的场合则是完成二者的估计,并在满足非负与和为一约束的前提下使得N尽可能小.

1.2 非负矩阵分解

非负矩阵分解是在分解后所有元素都非负的要求下将一个高维矩阵分解为两个低维矩阵的乘积,即

$ \mathit{\boldsymbol{V}} \approx \mathit{\boldsymbol{WH}} $ (3)

通常将n×m维矩阵V称为原矩阵,n×r维矩阵W称为基矩阵,r×m维矩阵H称为系数矩阵,其中r为特征数目,整个过程可以看作对原矩阵进行降维和特征提取.容易看出,非负矩阵分解的结果并不唯一,存在非负可逆矩阵D使得WDD-1H仍然是一个非负矩阵分解结果.因此需要定义一个目标函数来刻画分解后的乘积对原矩阵的逼近程度,目前最常用的是最小化二者之差的欧式距离平方,即Frobenius范数.Sajda则给出了非负矩阵分解的概率模型[18],目标为最小化损失函数,即

$ f\left( {\mathit{\boldsymbol{W}},\mathit{\boldsymbol{H}}} \right) = \frac{1}{2}\left\| {\mathit{\boldsymbol{V}} - \mathit{\boldsymbol{WH}}} \right\|_{\rm{F}}^2 = \sum\nolimits_{i,j} {{{\left[ {{\mathit{\boldsymbol{V}}_{i,j}} - {{\left( {\mathit{\boldsymbol{WH}}} \right)}_{i,j}}} \right]}^2}} $ (4)

最早使用传统梯度下降法和加性迭代规则求解该优化问题,随之而来的问题是步长难以控制以及迭代过程中非负性无法保障,在这方面Lee和Seung做出了开创性的工作,不仅提出了乘性迭代规则并且采用辅助函数的手段证明了其收敛[19],即

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{W}}_{i,k}} \leftarrow {\mathit{\boldsymbol{W}}_{i,k}}\frac{{{{\left( {\mathit{\boldsymbol{V}}{\mathit{\boldsymbol{H}}^{\rm{T}}}} \right)}_{i,k}}}}{{{{\left( {\mathit{\boldsymbol{WH}}{\mathit{\boldsymbol{H}}^{\rm{T}}}} \right)}_{i,k}}}}\\ {\mathit{\boldsymbol{H}}_{k,j}} \leftarrow {\mathit{\boldsymbol{H}}_{k,j}}\frac{{{{\left( {{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{V}}} \right)}_{k,j}}}}{{{{\left( {{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{WH}}} \right)}_{k,j}}}} \end{array} \right. $ (5)

相比加性迭代规则,乘性迭代规则具备很大的优势,可对步长进行精准控制,迭代过程中每个矩阵元素施以不同的步长而非对整个矩阵使用同一步长,从而保证矩阵元素不会出现负值,且很容易通过编码实现.

2 最大整体包容度约束的非负矩阵分解

传统的NMF算法在非负性限定之外,仅考虑分解后的整体逼近效果,而式(4)给出的目标函数非凸,因此存在多个局部最优解,导致将NMF用于具体场合时,往往已经达到最佳近似,分解效果却不尽如人意.应用NMF做高光谱图像像元分析时,需要加入额外反映端元某种特性的约束与Frobenius范数一起构成目标函数,并采用比例因子在分解的整体逼近和端元的近似程度之间取得平衡.约束NMF框架下的端元提取算法大多属于上述范畴,本章提出一个反映端元整体包容度的计算模型,在像元高度混合的局面下优先考虑端元的整体覆盖程度.

2.1 端元的最大整体包容度

图像中存在纯像元的理想状态下,所有的像元应该落于以端元为顶点的高维单形体内,处在光谱特征空间的一个子空间中.在高度混合无纯像元存在的情况下,将包含像元数量最多单形体的顶点作为端元的估计是一项合理选择.设端元集合E=[e1, e2, …, ep],p为端元个数,ein×1维端元向量,在p-1维线性空间E中构成单形体S.为方便计算,将0向量纳入端元集合,整体构成一个p维线性空间,物理上可以理解为纯黑背景也是一个端元,不影响最终的端元提取结果.在不降维的前提下计算该单形体的体积VS[20], 即

$ {V_{\rm{S}}} = \frac{1}{{p!}}\sqrt {\det \left( {{\mathit{\boldsymbol{E}}^{\rm{T}}}\mathit{\boldsymbol{E}}} \right)} $ (6)

将任一混合像元xt加入矩阵E,形成E+=[e1, e2, …, ep, xt],计算体积Vt

$ {V_t} = \frac{1}{{\left( {p + 1} \right)!}}\sqrt {\det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)} $ (7)

如果混合像元xt被包容在由端元构成的单形体内,则Vt=0,否则Vt>0.给出描述端元整体包容性的约束函数J(V)为

$ J\left( V \right) = \sum\nolimits_{t = 1}^m {{{\left( {{V_t}} \right)}^2}} $ (8)
2.2 算法与实现

经分析可知,J(V)值越小,表示端元集合对于所有像元的整体包容度越好,把J(V)作为新的约束项加入NMF框架中,给出最大整体包容度约束的非负矩阵分解(Maximum Overall Coverage Constraint Nonnegative Matrix Factorization, MOCC-NMF),式中J(V)是惩罚因子,用以控制对矩阵分解整体精确度和端元包容度之间的侧重,则

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{EC}}\;\;\;\mathit{\boldsymbol{E}} \ge 0,\mathit{\boldsymbol{C}} \ge 0\\ \left[ {\mathit{\boldsymbol{E}},\mathit{\boldsymbol{C}}} \right] = \begin{array}{*{20}{c}} {\arg \min }\\ {\mathit{\boldsymbol{E}},\mathit{\boldsymbol{C}}} \end{array}{f_{{\rm{MOCC - NMF}}}}\left( {\mathit{\boldsymbol{E}},\mathit{\boldsymbol{C}}} \right)\\ {f_{{\rm{MOCC - NMF}}}}\left( {\mathit{\boldsymbol{E}},\mathit{\boldsymbol{C}}} \right) = \frac{1}{2}\left\| {\mathit{\boldsymbol{X}} = \mathit{\boldsymbol{EC}}} \right\|_{\rm{F}}^2 + \mathit{\boldsymbol{\lambda J}}\left( \mathit{\boldsymbol{V}} \right) \\= \sum\nolimits_{i,j} {{{\left[ {{\mathit{\boldsymbol{X}}_{i,j}} - {{\left( {\mathit{\boldsymbol{EC}}} \right)}_{i,j}}} \right]}^2}} + \lambda \sum\nolimits_{t = 1}^m {{{\left( {{V_t}} \right)}^2}} \end{array} \right. $ (9)

求解MOCC-NMF的过程仍是基于传统的梯度下降思想,额外约束项的加入使得目标函数梯度计算发生变化,容易看出,fMOCC-NMF(E, C)中新加入的约束项与丰度矩阵C无关,易有

$ \left\{ \begin{array}{l} {\nabla _\mathit{\boldsymbol{C}}}{f_{{\rm{MOCC - NMF}}}}\left( {\mathit{\boldsymbol{E}},\mathit{\boldsymbol{C}}} \right) = {\mathit{\boldsymbol{E}}^{\rm{T}}}\mathit{\boldsymbol{EC}} - {\mathit{\boldsymbol{E}}^{\rm{T}}}\mathit{\boldsymbol{X}}\\ {\nabla _\mathit{\boldsymbol{E}}}{f_{{\rm{MOCC - NMF}}}}\left( {\mathit{\boldsymbol{E}},\mathit{\boldsymbol{C}}} \right) = \mathit{\boldsymbol{EC}}{\mathit{\boldsymbol{C}}^{\rm{T}}} - \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{C}}^{\rm{T}}} + \lambda \frac{{\partial J\left( V \right)}}{{\partial \mathit{\boldsymbol{E}}}} \end{array} \right. $ (10)

式中J(V)是一个连加式,其中每一项对E的偏导有

$ \frac{{\partial {{\left( {{V_t}} \right)}^2}}}{{\partial \mathit{\boldsymbol{E}}}} = \frac{1}{{{{\left[ {\left( {p + 1} \right)!} \right]}^2}}}\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial \mathit{\boldsymbol{E}}}} = \frac{1}{{{{\left[ {\left( {p + 1} \right)!} \right]}^2}}}\left[ {\begin{array}{*{20}{c}} {\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{1,1}}}}}&{\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{1,2}}}}}& \cdots &{\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{1,p}}}}}\\ {\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{2,1}}}}}&{\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{2,2}}}}}& \cdots &{\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{2,p}}}}}\\ \vdots&\vdots &{}& \vdots \\ {\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{n,1}}}}}&{\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{n,2}}}}}& \cdots &{\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{n,p}}}}} \end{array}} \right] $ (11)

注意到E+=[e1, e2, …, ep, xt],为方便书写,用符号ep+1表示xt,此处并不代表xt是一个端元.式(11)中矩阵的每一项有

$ \frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{i,j}}}} = \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right){\rm{tr}}\left[ {{{\left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}^{ - 1}}\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{i,j}}}}} \right] $ (12)

式中tr()为矩阵的迹.E+TE+是一个(p+1)×(p+1)维对称阵,且(E+TE+)i, j=〈ei, ej〉,将(E+TE+)-1用伴随矩阵的方式写出,Mei, ej〉为元素(E+TE+)i, j的代数余子式,有

$ {\rm{tr}}\left[ {{{\left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}^{ - 1}}\frac{{\partial \det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}{{\partial {\mathit{\boldsymbol{E}}_{i,j}}}}} \right] = \frac{2}{{\det \left( {\mathit{\boldsymbol{E}}_ + ^{\rm{T}}{\mathit{\boldsymbol{E}}_ + }} \right)}}\sum\nolimits_{l = 1}^{p + 1} {M\left\langle {{\mathit{\boldsymbol{e}}_j},{\mathit{\boldsymbol{e}}_l}} \right\rangle {{\left( {{\mathit{\boldsymbol{E}}_ + }} \right)}_{i,l}}} $ (13)

τ表示结合式(11)~(13)出现的常数项$ \frac{2}{{{{\left[{\left( {p + 1} \right)!} \right]}^2}}} $,给出适用于MOCC-NMF算法的乘性迭代规则,即

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{E}}_{i,k}} \leftarrow \frac{{{\mathit{\boldsymbol{E}}_{i,k}}}}{{{{\left( {\mathit{\boldsymbol{EC}}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right)}_{i,k}}}}\left[ {{{\left( {\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right)}_{i,k}} - \lambda \tau \sum\nolimits_{t = 1}^m {\sum\nolimits_{l = 1}^{p + 1} {M\left\langle {{\mathit{\boldsymbol{e}}_j},{\mathit{\boldsymbol{e}}_l}} \right\rangle {{\left( {{\mathit{\boldsymbol{E}}_ + }} \right)}_{i,l}}} } } \right]\\ {\mathit{\boldsymbol{C}}_{k,j}} \leftarrow {\mathit{\boldsymbol{C}}_{k,j}}\frac{{{{\left( {{\mathit{\boldsymbol{E}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right)}_{k,j}}}}{{{{\left( {{\mathit{\boldsymbol{E}}^{\rm{T}}}\mathit{\boldsymbol{EC}}} \right)}_{k,j}}}} \end{array} \right. $ (14)

设定迭代总次数或者误差阈值,按式(14)迭代求解即可完成MOCC-NMF.式(14)中求解Ei, k的过程中,有出现负值的可能性,届时可动态调整λ的大小,弱化约束项的作用,也可直接在该次迭代中将约束项置零,此次迭代则回归传统的NMF方法.实际上多次实验表明,带有此约束项的乘性迭代规则出现负值的概率并不大.

具体算法执行时还需要处理一些细节,在NMF中,低维特征数r和迭代初始值是很重要的两个参数,对于高光谱图像像元分析就是估计端元个数p和选取初始端元集E0.这里采用虚拟维度法(Virtual Dimensionality VD)[21]确定p,可以随机在图像中选择相关性不强的p个向量构成E0,也可用某种较简单的端元提取算法,比如N-FINDR等,生成E0.并通过对像元矩阵X和端元矩阵E底部分别加上m维和p维元素全为1的行向量来处理对于丰度的和为一约束.算法流程为:

1) 设定迭代次数或者误差阈值;

2) 估计端元个数p,生成初始端元集E0

3) $ \mathit{\boldsymbol{\overline X}} = \left[\begin{array}{l} \mathit{\boldsymbol{X}}\\ {\delta _1} \end{array} \right], {\mathit{\boldsymbol{\overline E}} _0} = \left[\begin{array}{l} {\mathit{\boldsymbol{E}}_0}\\ {\delta _1} \end{array} \right] $

4) 迭代开始;

5) $ {\mathit{\boldsymbol{E}}_{i, k}} \leftarrow \frac{{{\mathit{\boldsymbol{E}}_{i, k}}}}{{{{\left( {\mathit{\boldsymbol{EC}}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right)}_{i, k}}}}\left[{{{\left( {\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right)}_{i, k}}-\lambda \tau \sum\nolimits_{t = 1}^m {\sum\nolimits_{l = 1}^{p + 1} {M\left\langle {{\mathit{\boldsymbol{e}}_j}, {\mathit{\boldsymbol{e}}_l}} \right\rangle {{\left( {{\mathit{\boldsymbol{E}}_ + }} \right)}_{i, l}}} } } \right] $

6) $ {\mathit{\boldsymbol{C}}_{k, j}} \leftarrow {\mathit{\boldsymbol{C}}_{k, j}}\frac{{{{\left( {{\mathit{\boldsymbol{E}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right)}_{k, j}}}}{{{{\left( {{\mathit{\boldsymbol{E}}^{\rm{T}}}\mathit{\boldsymbol{EC}}} \right)}_{k, j}}}} $

7) 若迭代条件达成,输出EC,否则返回4).

3 实验与分析

分别在模拟图像和真实图像中对本文算法进行评测.参与结果比较的算法选用传统的N-FINDR和VCA端元提取算法,并采用完全约束的最小二乘法(Fully Constrained Least Squares Method, FCLS)[22]做丰度反演.无约束的NMF和MVC-NMF也参与实验比对.用光谱角距离(Spectral Angle Distance, SAD)和光谱信息散度(Spectral Information Divergence, SID)[23]衡量输出端元和参考光谱的相似度,用原丰度和估计丰度之间的均方根误差(Root Mean Square Error, RMSE)[24]评价丰度反演的效果.

3.1 模拟图像实验

从USGS数字光谱库中任选11条光谱作为纯像元生成模拟图像,光谱波段范围0.39~2.56 μm,光谱分辨率小于30 nm,波段数为420.模拟图像空间分辨率为105×105像素,端元分别处于图像的四角、中心点以及左右边界至中心的两条垂直分隔线上等比例的六个位置,每个纯像元各占5×5像素.混合像元处在图像被纯像元分割出的十个区域内,按线性光谱混合模型,丰度与相应端元距离成反比的规律进行混合,即离某端元越远,该端元在此像元中所占丰度越小,与现实中情况相吻合,且丰度满足和为一约束.然后去除所有包含丰度大于0.85项的像元,用相邻的混合像元替代,生成无纯像元的模拟图像.对于高度混合的图像,N-FINDR和VCA均没有唯一解,而传统NMF本就可能存在多个局部最优解,因此将N-FINDR、VCA和传统NMF算法多次运行(这里取200次),对输出做平均后参与比较.实验平台为i5-4950/12G/WIN7 64位,仿真环境为Matlab2015a.NMF框架下三个算法的终止条件均设为迭代总次数到达300时停止,MVC-NMF相关参数依据文献[13]设置.惩罚因子λ与实验图像以及约束项密切相关,故设定较为复杂.对于本文MOCC-NMF算法,给出规范化惩罚因子λ=λτ/det(E0TE0).根据多次实验结果并结合误差分析,得出结论λ取值处于区间[10-5, 10-4],并在λ≈3.784×10-5处分解的整体误差最小,因此对于本组实验数据MOCC-NMF算法的规范化惩罚因子λ围绕3.784×10-5设定.图 1给出了不同算法在相应指标下的表现.

图 1 不同算法在高度混合模拟图像下性能比较 Fig.1 Performance comparison of different algorithms using highly mixed simulated data

模拟图像中不存在纯像元,而N-FINDR和VCA本质上是端元选择算法,无法得到相对正确的输出,在相应指标下表现较差.传统的NMF算法直接用于混合像元分析,因为仅仅强调分解的整体相似性而不关注端元的任何其他特性,在没有好的初始值的帮助下,局部最优解距离实际最优结果有较大距离.MVC-NMF和本文算法从不同侧面考虑到了端元的某些特性,分别是特征空间中容纳全部像元前提下的最小单形体体积和端元对混合像元的包容程度,因此效果较其他算法好,而在高度混合的场景下,实验表明端元的整体包容性是更值得关注的特性之一.

在时间复杂度方面,尽管本文的MOCC-NMF算法在目标函数中加入的额外约束项涉及大量行列式计算,基于下述原因,算法的实际运行时间适中:首先,随着一个较好的初始端元集E0开始算法,E0由某种端元选择算法估计,从而本身具备较好的整体包容度,大多数引入混合像元形成新单形体体积为零,通过预判可避免大部分行列式的具体计算;其次,所加入额外的约束项可以使MOCC-NMF算法在更少的迭代步数中收敛.表 1给出了相同硬件配置下,不同算法在无噪声模拟图像上的运行时间.可以出看出,相比于传统的NMF以及加入约束项的同类算法,本文算法表现出不高的时间复杂度.此外,随着MOCC-NMF中主要计算部分,矩阵乘法和行列式计算并行化的实现,算法所耗时间将进一步降低.

表 1 不同算法在无噪声模拟图像上的运行时间比较 Tab.1 Comparison of different algorithms in elapse time using synthetic image without noise

为模拟图像加入零均值的高斯白噪声,用以评价算法在不同噪声环境下的鲁棒性,定义信噪比为SNR=10log10E[XTX]/E[nTn].图 2给出了算法在不同噪声情况中各个指标的表现.

图 2 算法在不同信噪比下的性能比较 Fig.2 Performance comparison under different SNR

信噪比的降低无疑会影响算法性能,MOCC-NMF在噪声干扰大的环境中会把噪声点当成普通像元试图包容在单形体中,从而给出错误的估计,通过对噪声的预处理可以有效避免此类问题.从图 2中可以看出,当信噪比大于25 dB,即接近实际高光谱图像信噪比上界时,本文算法性能开始优于其他算法.

3.2 真实图像实验

实验选用由机载可见光/红外光谱成像仪AVIRIS于1995年获取的美国内华达州南部Cuprite地区的高光谱遥感图像.Cuprite数据含224个波段,本次实验使用移除水吸收和低信噪比波段后的短波红外部分(波长1 990~2 479 nm),包含50个波段,空间分辨率为350×400,数据为经过大气校正的反射谱,由ENVI软件提供,光谱幅值均作了整数化处理.该地区地表多为裸露矿物,像元混合程度适中,且经过多次地面实际勘测具有较完备的地表信息,常被用于遥感地质光谱研究相关实验.结合ENVI软件给出的测试数据和VD估计的端元个数,去除噪声,合并重复地物,最终确定图像所反映区域包含9个端元,分别为沸石类(Zeolites)、黑色背景(Dark/background)、沙漠盆地(Playa)、方解石(Calcite)、明矾石(Alunite)、高岭石(Kaolinite)、白云母(Muscovite)、硅石(Silica)以及水铵长石(Buddingtonite).不同算法输出端元和参考光谱之间的光谱角距离(弧度)如表 2.

表 2 五种算法输出端元与参考光谱的光谱角距离(弧度)比较 Tab.2 Spectral angle distance (radian) between extracted endmembers and reference spectra by five algorithms

表中数据表明,NMF框架下的混合像元分析方法整体优于单纯的端元选择类算法,无论是否有纯像元存在于图像中,NMF类算法总给出相对准确的估计.不加约束的NMF算法因为局部最优解的问题无法得到满意的效果,即便是多次运行之后求平均仍是如此,而施加合适的约束之后效果明显得以提升.图 3给出了不同算法提取的光谱曲线与参考光谱曲线的对比.

图 3 不同算法提取的光谱曲线比较 Fig.3 Comparison of endmember spectra extracted by different algorithms

图 3可以看出,本文算法输出的端元光谱曲线与参考光谱最为接近,能够更好地反映端元地物在不同波段下的吸收特征.比如沸石类(Zeolites)端元在波长2.21 μm、2.34 μm左右处的吸收特性,白云母(Muscovite)端元在波长范围2.26~2.36 μm内的吸收特性,均在本文算法的输出中有较为明显的反映.

4 结论

基于NMF框架下的高光谱遥感图像混合像元分析思路,通过分析图像数据在特征空间中的凸面几何学性质,给出了反映端元整体包容性的计算方法,并据此提出了最大整体包容度约束的非负矩阵分解(MOCC-NMF)算法.考虑在NMF整体近似与非负性的同时,将包含像元数量最多作为新的约束,并通过惩罚因子调配约束与整体近似的关系.新目标函数中的约束项通过计算体积来表达端元的整体包容度,在计算过程中不对原始数据降维,因此不破坏数据在特征空间的表征,也不损害其物理意义.算法具体实现时,仍采用乘性迭代规则,充分利用其个体步长具备的优势,计算量方面,虽然体积的计算涉及行列式,但绝大多数被包容的像元约束项为零,执行效率不低,在模拟和真实图像上的相关实验证实了算法的有效性.

后续研究工作将在如下几个方面展开:1)寻找更适合反映高光谱遥感图像某些特性的新约束项;2)为NMF框架下整体近似与约束项之间关系的确定提供理论依据;3)低维特征数以及初始化是NMF重要的两个参数,将结合高光谱图像对二者展开研究;4)深入研究本文算法的数学特性、收敛性及更为适用的迭代规则;5)研究惩罚因子和算法以及数据的内在联系;6)实现相应的并行算法.

参考文献
[1] TONG Qing-xi, ZHANG Bing, ZHENG Lan-fen. Hyperspectral remote sensing:principle, technology and application[M]. Beijing: Higher Education Press, 2006.
童庆禧, 张兵, 郑兰芬. 高光谱遥感:原理、技术与应用[M]. 北京: 高等教育出版社, 2006.
[2] JIMENEZ L O, LANDGREBE D A. Supervised classification in high-dimensional space:geometrical, statistical, and asymptotical properties of multivariate data[J]. IEEE Transactions on Systems, Man, and Cybernetics, Part C:Applications and Reviews, 1998, 28(1): 39-54. DOI:10.1109/5326.661089
[3] NASCIMENTO J M P, BIOUCAS-DIAS J M. Hyperspectral unmixing based on mixtures of Dirichlet components[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(3): 863-878. DOI:10.1109/TGRS.2011.2163941
[4] KESHAVA N. A survey of spectralunmixing algorithms[J]. Lincoln Laboratory Journal, 2003: 55-78.
[5] PLAZA A, MARTÍN G, PLAZA J, et al. Recent developments in endmember extraction and spectral unmixing[M]. Optical Remote Sensing, 2011: 235-267.
[6] BOARDMAN J W, KRUSE F A, GREEN R O. Mapping target signatures via partialunmixing of AVIRIS data[C]. Summaries 5th JPL Airborne Earth Science Workshop 1995.1:23-26. http://www.oalib.com/references/18903744
[7] CHANG C I, PLAZA A. A fast iterative algorithm forimplementation of pixel purity index[J]. IEEE Geoscience and Remote Sensing Letters, 2006, 3(1): 63-67. DOI:10.1109/LGRS.2005.856701
[8] WINTER M E. A proof of the N-FINDR algorithm for the automated detection of endmembers in a hyperspectral image[C].Defense and Security. International Society for Optics and Photonics, 2004:31-41. http://spie.org/Publications/Proceedings/Paper/10.1117/12.542854
[9] ZHAO Chun-hui, GUO Yun-ting. An improved fast N-FINDR endmember extraction algorithm[J]. Acta Photonica Sinica, 2015, 44(10): 1022001.
赵春晖, 郭蕴霆. 一种改进的快速N-FINDR端元提取算法[J]. 光子学报, 2015, 44(10): 1022001.
[10] CRAIG M D. Minimum volume transforms for remotely sensed data[J]. IEEE Transactions on Geoscience and Remote Sensing, 1994, 32(3): 542-552. DOI:10.1109/36.297973
[11] NASCIMENTO J M P, DIAS J M B. Vertex component analysis:a fast algorithm tounmix hyperspectral data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(4): 898-910. DOI:10.1109/TGRS.2005.844293
[12] LIU Wei-xiang, ZHENG Nan-ning, YOU Qu-bo. Nonnegative matrix factorization and its applications in pattern recognition[J]. Chinese Science Bulletin, 2006, 51(1): 7-18. DOI:10.1007/s11434-005-1109-6
[13] MIAO Li-dan, QI Hai-rong. Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(3): 765-777. DOI:10.1109/TGRS.2006.888466
[14] LIU Jian-jun, WU Ze-bin, WEI Zhi-hui, et al. A fast algorithm for hyperspectral unmixing based on constrained nonnegative matrix factorization[J]. Acta Electronica Sinica, 2013, 41(3): 432-437.
刘建军, 吴泽彬, 韦志辉, 等. 基于约束非负矩阵分解的高光谱图像解混快速算法[J]. 电子学报, 2013, 41(3): 432-437.
[15] GENG Xiu-rui, JI Lu-yan, ZHAO Yong-chao, et al. A new endmember generation algorithm based on a geometric optimization model for hyperspectral images[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(4): 811-815. DOI:10.1109/LGRS.2012.2224635
[16] JI Lu-yan, GENG Xiu-rui, YU Kai, et al. A new non-negative matrix factorization method based on barycentric coordinates for endmember extraction in hyperspectral remote sensing[J]. International Journal of Remote Sensing, 2013, 34(19): 6577-6586. DOI:10.1080/01431161.2013.804223
[17] YANG Bin, LUO Wen-fei. Constrained NMF-based high-dimension adaptive particle swarm optimization algorithm for endmember extraction from a hyper spectral remote sensing image[J]. Journal of Remote Sensing, 2015, 19(2): 240-253.
杨斌, 罗文斐. 约束非负矩阵分解框架下高维自适应粒子群端元提取[J]. 遥感学报, 2015, 19(2): 240-253.
[18] SAJDA P, DU S. Recovery of constituent spectra using non-negative matrix factorization[C]. SPIE, 2003, 5207:321-331.
[19] LEE D D. Algorithms for nonnegative matrix factorization[J]. Advances in Neural Information Processing Systems, 2000, 13(6): 556-562.
[20] GENG Xiu-rui, ZHAO Yong-chao, WANG Fu-xiang, et al. A new volume formula for a simplex and its application to endmember extraction for hyperspectral image analysis[J]. International Journal of Remote Sensing, 2010, 31(4): 1027-1035. DOI:10.1080/01431160903154283
[21] CHANG C I, DU Q. Estimation of number of spectrally distinct signal sources in hyperspectral imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(3): 608-619. DOI:10.1109/TGRS.2003.819189
[22] 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.
[23] CHANG C, HEINZ D C. Constrainedsubpixel target detection for remotely sensed imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(3): 1144-1159. DOI:10.1109/36.843007
[24] PLAZA A, MARTINEZ P, PEREZ R, et al. A quantitative and comparative analysis of endmember extraction algorithms from hyperspectral data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(3): 650-663. DOI:10.1109/TGRS.2003.820314