一种应用于高光谱图像处理的非负矩阵分解方法

摘要

本发明公开一种基于稀疏性和相关性约束的非负矩阵分解方法(Non-negative Matrix Factorization，NMF)，并将该方法应用到高光谱遥感图像的混合像元分解的处理中。此方法最终将把给定的非负矩阵Vm×n分解成基矩阵Wm×r和系数矩阵Hr×n的乘积，即Vm×n≈Wm×rHr×n，首先选择非负矩阵V,随机初始化W和H，然后在目标函数中对系数矩阵H施加最小相关性约束，对基矩阵W和系数矩阵H施加稀疏性约束，然后按照迭代公式进行迭代至矩阵W和H收敛。

说明书

技术领域

本发明涉及一种高光谱图像处理方法，具体为一种基于稀疏性和相关性约束 的非负矩阵分解算法(Non-negative Matrix Factorization，NMF)，属于高光谱遥 感图像处理技术领域。

背景技术

遥感技术(Remote Sensing，RS)，指的是开始于二十世纪六十年代，从远距 离空间(航天)或者外太空空间(航空)对地球表面进行的综合性的观测技术。 遥感，即从远距离或外太空间接观测监控，不接触目标，从光学角度获取对目标、 现象以及区域的相关信息，从而进行数据的融合分析以及推断，最终达到获取所 需的目标信息的一种手段、技术和科学。高光谱图像遥感技术(Hyper-spectral  Graphic Remote Sensing)，是一种立于高光谱技术的技术基础上的遥感探测技术， 是一种融合新型光谱探测技术、微信号探测技术、光学精密机械、信号高速处理 技术、计算机处理信息技术在一起的先进性、综合性、科学性技术。该技术同时 由于影响到地球学、环境保护学、野生生物学、信息技术学、空间地理学等众多 科学领域的关键技术，其技术的进步发展备受国内外学者的强烈瞩目，目前已广 泛应用于植被绿化、土壤分析、精细农业管理、大气环境监测、水环境监测、勘 探矿产资源分布等方面，充分展示出了高光谱遥感技术的潜力与优势。

在全球遥感界中通常会这样认为：多光谱遥感(Multi-spectral Remote Sensing) 指的是光谱辨析率在1～9.9×10^-1λ的光谱范围内的遥感，在这种范围内的遥感传 感器内只有分布在可见光区和近中红外光光谱区的很少的波段数；高光谱遥感 (Hyper-spectral Remote Sensing)指的是光谱辨析率在1～9.9×10^-2λ的光谱范围 内的遥感技术；而超高光谱遥感(Ultra-spectral Remote Sensing)则指的是经过更高 技术的进步后，光谱辨析率达到1～9.9×10^-1λ的光谱范围内的遥感技术。遥感技 术，其发展历程是在全色(即黑白色)摄影与彩色影像阶段之后，多光谱遥感迅 速发展于二十世纪中后半叶时期，取得了较大的进步，并已运用于环境与资源探 测领域，但其分辨率仍然处于10^-1λ数量级光谱的波长范围内，只有相对较少的 采样点，随着科技的发展，人们对地球资源与环境认识越来越深入，其运用精度 越来越不能满足需求，对分辨率更深入的发展要求更加迫切，主要体现在空间分 辨率和光谱分辨率的深入与提高。图片中在空间分辨率指的是在遥感图片中相邻 的两个地形地物之间能够被区分出的最小长度，也就是通常所说的图片的清楚程 度，可用于对图片数据的直观解释。而图片中在光谱分辨率指的是光谱检测中的 光谱反射出的能够区分出的最小光谱间隔的波长，更确切的讲，就是检测区分光 谱的性能。相对于空间分辨率的直观观察，光谱观测技术能够更加直接和有效的 反应地物的结构与性质，尤其是在远洋探测监测、动植物研究分类、农业精细化 规划、勘探矿产调查和现代化军事运用方面具有更好的效果和更快的速率。随着 科技的日新月异，航天航空工业的飞速发展，遥感平台和光学传感器的不断更新 进步升级，提高遥感光谱分辨率的已经处在迅猛的发展势头中，高光谱遥感技术 的发展已经成为当今遥感技术研究学者们的研究热点。由于种类性质不一样的物 质在一定的波长的光谱照射下，具有各自不同比例的吸收和反射的特性，通过对 比其反射光谱(或吸收)光谱的之间的差异，可以导出其物质的组成成分和物理 结构上之间的差异。如何才能从数据复杂的图像中提取出物质的特有属性(特征) 并且清晰了解与周边物质及在整体形态中的关系，成为研究中主要解决的问题。

高光谱遥感影像(Hyper-spectral Remote Sensing Images)，是指在电磁波谱 中的中远红外光谱、近红外光谱、可见光光谱和紫外光谱区域内，在光谱成像仪 的作用下，得到的很多光谱分布连续且光谱区域范围很小的影像信息数据(如图 1所示)。随着成像光谱仪的迅速发展与进步，所能获取到的高光谱遥感影像相 比于传统的二维成像技术的遥感图像，其是具有几十以至上百个波段叠加而成的 光谱图，其中的每一个基本的像元结构都是从几十以至上百个连续波段通道中所 获取的得到的光谱图像，其正对应的光谱实物的反射特性，最后都会得到一条较 为完整的光谱曲线。高光谱遥感图像已不再是二维成像技术，将独特的光谱理论 与遥感成像技术有机的结合起来，形成的连续的地物光谱曲线，使得利用高光谱 技术能够成功反演地形地物的细节。

高光谱遥感图像技术是一种新型的综合性对地观测技术，该技术的发展已在 众多领域得到广泛的运用。高光谱图像由于其光谱分辨率高，在成像中能够记录 几十以至几百个波段的光谱信息，能够很好的识别分析地形地物的种类与性质， 其良好的效果与广泛的用途使得高光谱技术得到了众多学者们的格外关注，引起 了对高光谱图像技术的研究热潮。

相对于纯像元而言，混合像元是高光谱遥感图像中的每个数据单元是由不同 地物的光谱信息组成的混合信息。在现实中，由于地形地物分布不均、交综复杂 以及成像光谱仪的分辨率有限，导致高光谱遥感图像中一定会大量的存在着不均 匀像元(即混合像元)，尤其是在地物种类繁多且分布复杂的地形地域，若是将 图片中像元看成是纯净地物，那么会造成较大的误差和分析的精度下降。这使得 在处理分析过程中对混合像元的研究成为高光谱遥感图像的主流研究方向，成为 一个新的研究热点。混合像元的存在阻碍了对高光谱图像的研究利用，对混合像 元的分解成为了高光谱图像技术的基础。

1999年，D.D.Lee与H.S.Seung在美国著名科学杂志《Nature》上发表了一 篇极其重要的文章，其创造性地提出了一种带有非负性约束的矩阵分解的方法： 非负矩阵分解算法(Non-negative Matrix Factorization，NMF)。非负矩阵分解， 是一种在非负约束下的矩阵分解算法，由于其算法的便捷性和对广泛物理意义的 表征性使得非负矩阵分解算法得到广泛的运用。该算法是在非负性的约束下，将 原矩阵分解成一个基矩阵和系数矩阵的乘积，其中基矩阵和系数矩阵中所有元素 满足非负的约束。光谱信息的非负性也适合于非负矩阵的要求，分解后的结果是 端元向量矩阵与其丰度矩阵的乘积，在物理意义的吻合使得非负矩阵分解算法在 混合像元分解的研究中得到了广泛的运用。但是非负矩阵分解算法具有迭代收敛 慢和容易陷入局部最优解的局限性。

2002年P.O.Hoyer在IEEE(Institute of Electrical and Electronics Engineers) 上发表了重要文章，提出了一种将稀疏编码与NMF算法相融合的方法，该方法 是在NMF算法加入了一种稀疏度约束。W.Liu等人在此基础上提出了稀疏NMF (Sparse NMF，SNMF)算法。SNMF算法在NMF的目标函数中加入了关于表征 丰度稀疏性的约束项，而对于自然信号，在科学与工程领域中而言，稀疏性是一 种比较常见的属性。

加权NMF(Weighted NMF，WNMF)算法是众多的改进算法中的一类。它 是在NMF目标函数中加入了对角矩阵U，该矩阵中对角线上的元素值大小表示 着各个训练样本向量的重要程度，从而对各类样本的重要程度进行了区分，这就 会有以利于防止出现相似的训练样本向量会产生冗余信息。加权式的改进措施可 以适用于任何一个优秀的方法中。

在NMF分解中，产生的基矩阵W通常具有显示训练样本特征的物理意义， 而这特征也就是训练样本的局部性，进而局部NMF(Local NMF，LNMF)算法 就应运而生，该算法是在目标函数里加入了一定的约束因子，如此确保分解结果 尽可能的正交化，从而使得分解产生的结果能够更多的体现出训练样本的局部特 征。

多年来，在基本NMF算法的基础上，很多改进的NMF算法被提出，SNMF、WNMF、 LNMF等，不同的改进算法有着不同的运用方向，经实际运用证明，各改进算法 均具有更优异的不同的属性和更加令人满意的收敛速度，对于NMF算法的实际运 用取得了很大的增强效果。

发明内容

发明目的：在高光谱遥感图像的研究中，由于地形地物的复杂性，即使较高 的分辨率也不能使得图像的像元为纯净的像元，而是由像元空间区域内的各种端 元的混合反射光谱，即混合像元。混合像元的存在阻碍了对高光谱图像的研究利 用，非负矩阵分解是一种在非负约束下的有效的矩阵分解算法。为了利于高光谱 遥感图像中混合像元分解，使得分解所得的基矩阵和系数矩阵能够得到期望的高 稀疏性和低相关性的要求，对非负矩阵分解过程的目标函数中的基矩阵和系数矩 阵施加稀疏性约束并对系数矩阵施加相关性约束。

技术方案：一种应用于高光谱图像处理的非负矩阵分解方法，该方法考虑矩 阵的稀疏特性和相关性，并一定程度上影响了混合像元分解的精度。具体如下： 此方法最终将把给定的非负矩阵V_m×n分解成基矩阵W_m×r和系数矩阵H_r×n的乘积， 即V_m×n≈W_m×rH_r×n，首先选择非负矩阵V,随机初始化W和H，然后在目标函数 中对系数矩阵H施加最小相关性约束，对基矩阵W和系数矩阵H施加稀疏性约 束，然后按照迭代公式进行迭代至矩阵W和H收敛。

具体包括如下步骤：

步骤1：选择非负矩阵V，随机初始化W和H；

步骤2：通过对目标函数F添加了对系数矩阵H进行了最小相关性约束和对 基矩阵W和系数矩阵H都施加了2-范数约束，得到式子(5)，利用基本NMF 算法中的下降法给出一种更新规则，得到新的迭代公式(6)，按式子(6)进行 迭代；

$F(W,H)=\frac{1}{2}{{||V-WH||}_F}^2-\frac{1}{2}{\alpha }_W{{||W||}_F}^2-\frac{1}{2}{\alpha }_H{{||H||}_F}^2+{\beta }_{H}R\left(H\right)---\left(5\right)$

$\left(\begin{array}{c}W\leftarrow W./\left({\mathrm{WHH}}^T\right)·*({\mathrm{VH}}^T+{\alpha }_{W}W)\\ H\leftarrow H./\left(W^{T}WV\right)·*\left[W^{T}V-{\alpha }_{H}H-{\beta }_H\left(\frac{1}{\left(h_k{h_k}^T\right)H}-{\left({\mathrm{HH}}^T\right)}^{-T}H\right)\right]\end{array}\right)---\left(6\right)$

式中${{||W||}_F}^2=\underset{i,k}{\Sigma }{W_{ik}}^2,{{||H||}_F}^2=\underset{k,j}{\Sigma }{H_{kj}}^2;$α_W和α_H为保证迭代收敛的正则化参数， β_H是用来平衡NMF的精度和基矩阵W以及系数矩阵H特性的参数。

步骤3：在每一步迭代后将W和H的负元素置零，并对W按列进行归一化 处理；

步骤4：循环进行步骤二～三直至收敛，此时得到基矩阵W和系数矩阵H。

为了更好的理解本发明所涉及的技术和方法，在此对本发明涉及的理论进行 介绍。

1、基本NMF算法

假设给定原矩阵V_m×n满足非负，则NMF算法的目的是求解出两个非负低秩 子矩阵W_m×r和H_r×n，从而使得

V_m×n≈W_m×rH_r×n      (7)

其中W_m×r称为基矩阵或特征矩阵，H_r×n称为系数矩阵。一般情况下 r＜＜min(m,n)，从而达到降低数据存储维数的目的，同时使得子矩阵W_m×r和H_r×n具有更加明确的解释意义。

NMF算法是一种低秩逼近算法，通常其目标函数可以是多种多样的， D.D.Lee与H.S.Seung在文章中主要对以下两种目标函数进行了研究推导：

一是运用原矩阵V和基矩阵W、系数矩阵H之间的欧氏距离(Euclidian  Distance)，将其平方作为目标函数，推导出迭代公式，然后通过逐步的迭代，尽 量使得欧氏距离最小化从而达到最佳逼近的目的，其目标函数为：

$F(W,H)=\frac{1}{2}{{||V-WH||}_F}^2=\frac{1}{2}\underset{i,j}{\Sigma }{(V_{ij}-\underset{k=1}{\overset{r}{\Sigma }}{W}_{ik}{H}_{kj})}^2---\left(8\right)$

当F(W,H)＝0时，V＝WH。目标函数式(8)关于W和H分别都是凸的， 但对W和H联合起来时并非是凸的，即同时关于W和H的时候目标函数是非凸 的，若在全局范围内求解最小值则是非常困难的，所以NMF求出的结果只能是 局部最优解。利用梯度下降法给出了一种更新规则，首先根据目标函数式(8) 对W和H分别求梯度得到：

$\{\begin{array}{c}\frac{\partial F}{\partial W_{ik}}=-\underset{i,k}{\Sigma }{\left(V^{T}H\right)}_{ik}+\underset{i,k}{\Sigma }{\left({\mathrm{WHH}}^T\right)}_{ik}\\ \frac{\partial F}{\partial H_{kj}}=-\underset{k,j}{\Sigma }{\left(W^{T}V\right)}_{ik}+\underset{k,j}{\Sigma }{\left(W^{T}WH\right)}_{kj}\end{array}---\left(9\right)$

再通过逐步迭代，其迭代规则为：

$\left(\begin{array}{c}{W}_{ik}\leftarrow W_{ik}-{\eta }_W\frac{\partial F}{\partial W_{ik}}\\ H_{kj}\leftarrow H_{kj}-{\eta }_H\frac{\partial F}{\partial H_{kj}}\end{array}\right)---\left(10\right)$

上式中当学习速率η_W、η_H按下式进行取值时：

$\{\begin{array}{c}{\eta }_W=\frac{W_{ik}}{{\left({\mathrm{WHH}}^T\right)}_{ik}}\\ {\eta }_H=\frac{Hkj}{{\left(W^{T}WH\right)}_{kj}}\end{array}---\left(11\right)$

将式(9)和式(11)分别代入式(10)，从而可以得到一个乘性迭代的规则：

$\{\begin{array}{c}{W}_{ik}\leftarrow W_{ik}\frac{{\left({\mathrm{VH}}^T\right)}_{ik}}{{\left({\mathrm{WHH}}^T\right)}_{ik}}\\ H_{kj}\leftarrow H_{kj}\frac{{\left(W^{\text{T}}V\right)}_{kj}}{{\left(W^{T}WH\right)}_{kj}}\end{array}---\left(12\right)$

将其写成矩阵的形式为：

$\{\begin{array}{c}W\leftarrow W.*\left({\mathrm{VH}}^T\right)./\left({\mathrm{WHH}}^T\right)\\ H\leftarrow H.*\left(W^{T}V\right)./\left(W^{T}WH\right)\end{array}---\left(13\right)$

二是运用V和WH之间的散度(KL-散度)，也就是概率意义上的距离作为目 标函数来达到最佳逼近：

$F(W,H)=\underset{i,j}{\Sigma }\left(V_{ij}\log \frac{V_{ij}}{{\left(WH\right)}_{ij}}-V_{ij}+{\left(WH\right)}_{ij}\right)---\left(14\right)$

该解法如同欧氏距离，可得到以下乘性迭代规则

$\left(\begin{array}{c}{W}_{ik}\leftarrow W_{ik}\frac{\underset{j}{\Sigma }{H}_{kj}{V}_{ij}/{\left(WH\right)}_{ij}}{\underset{j}{\Sigma }{H}_{kj}}\\ H_{kj}\leftarrow H_{kj}\frac{\underset{j}{\Sigma }{W}_{ik}{V}_{ij}/{\left(WH\right)}_{ij}}{\underset{i}{\Sigma }{W}_{ik}}\end{array}\right)---\left(15\right)$

NMF算法是一种低秩迭代算法，其实现的基础是目标函数的收敛性，下面 进行收敛性证明：

首先引入相关的定义，假定目标函数是F，则设G(h,h_t)为F(h)的辅助函数， 其满足条件：G(h,h_t)≥F(h)，当h＝h_t时相等，即G(h_t,h_t)＝F(h_t)，其中h为一 变量。G(h,h_t)与F(h)的关系如图2所示。

引理1：假定G为辅助函数，则在以下的迭代规则，F是非增的：

$h_{t+1}=\arg \underset{h}{\min }G(h,h_t)---\left(16\right)$

从而可以推出F(h_t+1)≤G(h_t+1,h_t)≤G(h_t,h_t)＝F(h_t)，因为当F(h)在h＝h_t为局部 最小值时，F(h_t+1)＝F(h_t)，则▽F(h_t)＝0，对于F(h)，根据迭代可以推出 $h_{\min }=\arg \underset{h}{\min F}\left(h\right),$即：

F(h_min)≤…≤F(h_t+1)≤F(h_t)≤…≤F(h₁)≤F(h₀)        (17)

对于利用目标函数对矩阵W和H的迭代过程是否收敛，我们只要能够构造 出上面所说的辅助函数，即可证明迭代过程收敛，下面以矩阵H为例，W的证 明与H是一样的，不再赘述。

引理2：假定K(h_t)为一个对角矩阵，K_ab(h_t)＝δ(W^TWh_t)_a/h_ta，其中δ函数 是一个冲激函数，则

$G(h,h_t)=F\left(h_t\right)+(h-h_t)▿F\left(h_t\right)+\frac{1}{2}{(h-h_t)}^{T}K\left(h_t\right)(h-h_t)---\left(18\right)$

其中F(h)为目标函数对矩阵H的迭代，即要证 明该G(h,h_t)函数为F(h)的构造函数只需证明G(h,h_t)≥F(h)与G(h_t,h_t)＝F(h_t) 即可。显然，把h＝h_t带入式(18)中G(h_t,h_t)＝F(h_t)，下面证明G(h,h_t)≥F(h)：

对F(h)在h＝h_t处进行泰勒公式展开即

$F\left(h\right)=F\left(h_t\right)+(h-h_t)▿F\left(h_t\right)+\frac{1}{2}{(h-h_t)}^T\left(W^{T}W\right)(h-h_t)---\left(19\right)$

要证明G(h,h_t)≥F(h)，只需证明G(h,h_t)-F(h)≥0，即

(h-h_t)^T(K(h_t)-W^TW)(h-h_t)≥0      (20)

为了证明其非负性，构造矩阵M_ij(h_t)＝h_ti(K(h_t)-W^TW)h_tj，只要证明M_ij(h_t) 是半正定的即可证明式(20)，即只要证明v^TMv≥0，即

$\left(\begin{array}{c}{v}^{T}Mv=\underset{ij}{\Sigma }{v}_i{M}_{ij}{v}_j\\ =\underset{ij}{\Sigma }\left[h_{ti}{\left(W^{T}W\right)}_{ij}{h}_{tj}{v_i}^2-v_i{h}_{ti}{\left(W^{T}W\right)}_{ij}{h}_{tj}{v}_j\right]\\ =\underset{ij}{\Sigma }{\left(W^{T}W\right)}_{ij}{h}_{ti}{h}_{tj}\left(\frac{1}{2}\underset{ij}{\Sigma }{v_i}^2+\frac{1}{2}{v_j}^2-v_i{v}_j\right)\\ =\frac{1}{2}\underset{ij}{\Sigma }\left(W^{T}W\right)h_{ti}{h}_{tj}{\left(v_i-v_j\right)}^2\\ \ge 0\end{array}\right)$得证。

所以，NMF算法的主要步骤和要求为：

(1)初始化W和H；

(2)选择学习速率和参数，保证迭代过程始终非负；

(3)按照迭代公式进行迭代；

(4)将W和H强制非负，且将H进行归一化处理；

(5)反复进行3、4至收敛。

2、稀疏NMF算法

2002年P.O.Hoyer在IEEE(Institute of Electrical and Electronics Engineers) 上发表了重要文章，提出了一种将稀疏编码与NMF算法相融合的方法，该方法 是在NMF算法加入了一种稀疏度约束。W.Liu等人在此基础上提出了稀疏NMF (Sparse NMF，SNMF)算法。所谓稀疏性，指的是矩阵中的元素数据分布的一种 特性，具有稀疏性的矩阵叫做稀疏矩阵，是指其矩阵中的元素中存在大量的零元 素且占有比例很高的矩阵。对于自然信号，在科学与工程领域中而言，稀疏性是 一种比较常见的属性。

SNMF算法是在NMF的目标函数中加入了关于表征丰度稀疏性的约束项， 从而使得目标函数变为

$F(W,H)=\frac{1}{2}{||V-WH||}^2+\lambda \underset{i,j}{\Sigma }{H}_{ij}---\left(21\right)$

约束项是以系数矩阵H中所有元素的代数和来表征对其稀疏度的约 束，并在迭代中使该约束项的数值达到最小，以使得结果尽可能的接近原矩阵， λ称为稀疏因子。

3、稀疏性约束

一般情况下，高光谱遥感图像中的各个像元中的端元分布不可能在整张图像 中分布，一般会比较集中分布于一片或若干片较为集中的区域。所以，在对高光 谱遥感图像非负分解的过程中施加稀疏性的约束，是符合实际情况，并且有利于 图像处理的后续工作的进行。

稀疏度，在矩阵论中对向量来说有一种利用1-范数与2-范数之间的关系作 为度量的方式，其具体方式为：

$S\left(x\right)=\frac{\sqrt{n}-\left(\Sigma \right|x_i\left|\right)/\sqrt{{{\mathrm{\Sigma x}}_i}^2}}{\sqrt{n}-1}---\left(22\right)$

其中n为向量x的维度，||x||₁＝Σ|x_i|为向量x的1-范数，为向 量x的2-范数。在对范数的学习中有这样的关系式：

$0\le \frac{1}{\sqrt{n}}{||x||}_1\le {||x||}_2\le {||x||}_1---\left(23\right)$

所以可以得到0≤S(x)≤1，当且仅当向量x只有一个元素时S(x)＝1，且当 向量x中所有元素相等且不为零时S(x)＝0。

可以看出，我们想得到较大稀疏度的条件是使得1-范数与2-范数之比较大， 即该式的值较小，因为在NMF算法的迭代过程中，对迭代的结果都 要作归一化处理，即1-范数的值等于1(||x||₁＝1)，所以我们要最大化2-范数， 可以使得分解后的结果具有较大的稀疏度。

在NMF算法中，在添加稀疏约束时不能仅仅是简单的从向量变成矩阵的2- 范数，矩阵的2-范数的定义为若仅仅是最大化基矩阵W的 2-范数，则有可能导致产生分解结果中W中有一个列向量的2-范数很大，但其他 的却很小，从而导致结果不能达到期望的稀疏性的要求，所以在施加约束是考虑 将最大化基矩阵W所有列向量的2-范数的平方之和作为目标，即

${{||W||}_F}^2=\underset{i,j}{\Sigma }{w_{ij}}^2---\left(24\right)$

同理对系数矩阵H施加为：

${{||H||}_F}^2=\underset{i,j}{\Sigma }{h_{ij}}^2---\left(25\right)$

4、相关性约束

相关性，指的是度量两个或两个以上的具备相互之间有所关联的性质的元素 变量之间的相互关联的大小指数。在高光谱图像中，各个像元的特征曲线之前的 相关程度主要存在于地物之间的差异，高光谱由于其光谱分辨率比较高，因而其 相关性应该比较小，但是由于地物的复杂性，混合像元在去光谱特征具有一定的 相似性，表现为具有一点的相关性，而各个端元之间的相关性是很小的。

在数理统计理论中，相关系数是指表现两个或两个以上的信号之间的相关性 的大小。相关系数越大表示相关性就越大，且在不相关的一组信号数据间所表现 出的相关系数的绝对值必定会小于它们之间的混合信号数据之间所表现出的相 关系数的绝对值。在NMF算法中，我们希望得到的结果是具有较小的相关性， 若在分离时得到的信号之间的相关系数的绝对值越小就表明分离出来的端元之 间通常都会包含一定量的相似的物质，在光谱特征上就变现为具有一定的相关 性。这种相关性通常会与像元之间的平滑性较为接近，表现为混合像元的光谱号 间的相关性越小，因此可以把分离信号相关系数最小作为对目标函数的迭代过程 中的一个约束条件。在信息学中，相关系数的定义为：

$R\left(H\right)=\frac{1}{2}\left[\underset{i=1}{\overset{n}{\Sigma }}log{\left({\mathrm{HH}}^T\right)}_{ii}-log\left|{\mathrm{HH}}^T\right|\right]---\left(26\right)$

当i≠j，若＜h_i,h_j＞＝0，则R(H)取最小值。

附图说明

图1为高光谱遥感图像示意图；

图2为辅助函数示意图；

图3为USGS数据库图片；

图4为灰度图像；

图5为南京禄口机场附近区域的实验结果图；

图6为南京青龙山附近区域的实验结果；

图7为南京江心洲附近区域的实验结果；

图8为南京栖霞区附近区域的实验结果；

图9为试验结果信噪比与维度r的关系；

图10为实验迭代次数与SNR的关系图。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本 发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发 明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

应用于高光谱图像处理的非负矩阵分解方法，包括如下步骤：

步骤1：选择非负矩阵V，随机初始化W和H；

步骤2：按式子(6)进行迭代；

步骤3：在每一步迭代后将W和H的负元素置零，并对W按列进行归一化 处理；

步骤4：循环进行步骤二～三直至收敛，此时得到基矩阵W和系数矩阵H。

仿真实验结果分析

1.实验图像

通过仿真实验对算法的性能进行分析和评价。为了验证基于约束的NMF算 法的有效性，本文对针对某一标准高光谱遥感影像进行试验。所使用的遥感影像 采用USGS数据库图像。

以下4张高光谱遥感图像是从USGS官网上下载的南京地区周围的4个区域 的图像。图像格式大小为420×314像素，因为高光谱遥感图像中的像元都混合的 端元组成的复杂光谱曲线，需要经过预处理，本实验采用USGS官网中可供下载 的已处理过的图片，可直接用于下一步处理。实际图片如图(3)所示。

USGS波普数据库其中主要是分析立体三维的成像光谱遥感数据，其包含了 多达498条，其共有有针对性地反应各种自然矿物质(即矿物矿产)资源的444 个植物或其他示例样品的共计498个反射波谱段数据。数据库的光谱波段范围为 0.2～3.0μm，波谱经过进一步的校正后，其实际波谱范围为0.4～2.5μm，其分 辨率在4nm(其波谱在可见光波段0.2～0.8μm)和10nm(其波谱在近红外波段 0.8～2.35μm)。

2、实验预处理和参数选择

下载的图片格式为彩色图片，大小为420×314×3，对于计算起来有些不方便， 所以首先转化成灰阶图片，其格式大小为420×314，像素值大小为0～1，为double 型。转化成的图像如下图4所示。

对该4张图片使用本算法进行试验，对参数的选择主要有降维后的基矩阵W 的列数和系数矩阵H的行数r、迭代算法的学习速率η_W和η_H以及世界约束条件 时的正则化参数α_W和α_H。经过多次试验对比，为了使得试验效率更高和结果更 好，本实验的参数选择为r＝140，学习速率取式(11)，正则化参数取α_W＝0.3， α_H＝0.5，β_H＝0.2，迭代次数为1000。

3、实验结果与分析

实验结果效果图如图5～8所示。其中V是原图像，WH是分解后重构的图像， 为了作以对比V-WH是指的重构前后的误差值，下面的W和H是分解得到的基 矩阵和系数矩阵。在这里我们主要关心的是误差值与分解的结果的优良性，但是 从图中我们无法辨别。为此，在试验中除了采用矩阵的稀疏度作为对比的标准外， 还将引入一个表征数据准确度的标准：信噪比(Signal-to-Noise Ratio，SNR)。 信噪比，又名“讯噪比”，指的是在电子领域中，某个设备或系统中的原信号与 在该设备或系统运行过程中所产生的噪声的比例大小。在图像处理的领域中，信 噪比有如下定义(带入本实验的变量)：

$SNR=10{\log }_{10}\left[\frac{{||V||}^2}{{||WH-V||}^2}\right]---\left(27\right)$

其中V为原图像的数据矩阵，WH为分解后重构的图像的数据矩阵，信噪比 SNR的数据量的单位采用分贝(dB)。在试验中我们将本文算法与NMF和SNMF 两种算法进行对比，结果如下：

表.1 南京禄口机场附近区域

算法 H的稀疏度 信噪比 NMF 0.2648 23.6317 SNMF 0.6184 25.6542 本文算法 0.6257 31.5333

表.2 南京青龙山附近区域

算法 H的稀疏度 信噪比 NMF 0.2896 23.4918 SNMF 0.6254 25.2826 本文算法 0.6458 31.1413

表.3 南京江心洲附近区域

算法 H的稀疏度 信噪比 NMF 0.3152 23.0784 SNMF 0.7185 25.3059 本文算法 0.7148 31.1424

表.4 南京栖霞区附近区域

算法 H的稀疏度 信噪比 NMF 0.2800 23.6474 SNMF 0.6605 25.2716 本文算法 0.6572 31.0521

通过以上实验结果可以明显发现本文算法会优于NMF和SNMF算法，与 NMF算法相比，在系数矩阵H的稀疏度上，本文算法占有很大优势，信噪比方 面有一定提升，与SNMF算法相比，在系数矩阵H的稀疏度上，与本文算法比 较接近，但是在信噪比方面，本文则会优秀很多。

在实验过程中，发现对于降维后的维数r并没有较多分析，仅仅只有 r＜＜min(m,n)的限定，而在试验中进行了对比与分析，参数选择如前所述，为 了达到实验效果，本次试验采用固定迭代次数，取1000以确保实验效果。实验 结果如图9所示。

从结果中可以看出当r＜50时，信噪比上升较快，此时的由于低维度r的取 值过低，导致原图像信息在分解的结果中信息丢失率比较严重，此时的r不适合 于试验取值；当100＜r＜200时，信噪比随低维度r的变化较为平稳，可以用此 段的r的值作为试验取值；当r＞200时，信噪比已基本不随低维度r的变化而有 较大变化了，且对于降维的效果也影响较大。试验中采用r＝140是比较合适的， 有一个不错的信噪比，即分离后重构的图像的效果比较好，同时，140相对于420 和314都小很多，说明了降维效果也比较好。总体能够达到我们的要求。

上述试验中迭代次数选为1000，使得迭代效果比较好，下面我们针对迭代 过程中的收敛条件改为$\frac{trace\left((WH-V){(WH-V)}^T\right)}{trace\left({\mathrm{VV}}^T\right)}\le 0.001,$函数trace()为矩阵对 角元素之和，得到的实验结果为：

表.5 试验收敛时信噪比与迭代次数

实验图片 实验结果信噪比 实验迭代次数 a 30.0018 454 b 30.0021 507 c 30.0052 442 d 30.0006 493

从实验结果表中可以发现，实验信噪比非常接近，究其原因，是因为我们设 定的收敛条件决定的，对比信噪比公式(26)与收敛条件可以发现信噪比是原矩 阵与误差矩阵(即WH-V)的行列式的比值加以运算的值，而我们收敛条件是 误差矩阵(即WH-V)与原矩阵的与其转置矩阵的乘积的迹的比值。由矩阵论 知识我们可以知道某方阵的行列式与其方阵乘以其转置的迹是成一定的比例。由 收敛条件中的系数值0.001，可以推导出信噪比 $SNR=10{\log }_{20}\left[\frac{{||V||}^2}{{||WH-V||}^2}\right]=10{\log }_{10}{\left[\frac{trace\left(\left(WH-V\right){\left(WH-V\right)}^T\right)}{trace\left({\mathrm{VV}}^T\right)}\right]}^{-1},$当 $\frac{trace\left((WH-V){(WH-V)}^T\right)}{trace\left({\mathrm{VV}}^T\right)}\le 0.001$时迭代收敛，即 时结束迭代，带入SNR计算式中，得 SNR＝10log₁₀(0.001)^-1＝30，所以我们的试验结果中信噪比都在30左右。实验 结果中迭代次数在450～500次左右，前面的实验中取迭代次数1000则显得有些 多余，但并不影响结果的正确性。为了研究迭代次数与信噪比的具体关系，下面 进行深入的实验，实验结果如下图10所示。

从图10可以看出迭代次数在400以下时，信噪比上升的速率很快，也就意 味着在这期间由于迭代次数不够导致分解的结果不太理想，重构的图片与原图片 差距较大，不能作为实验结果，当达到400至500次时信噪比上升的速率明显放 缓，迭代次数的上升对于信噪比的影响不是很大，当迭代次数超过500次时，基 本信噪比就不再变化了，已经接近实验比较好的结果了，此时再上升迭代的次数 已经没有什么收益了，也不适宜实验的效率。按照前面实验中得出450～500次 的迭代次数是比较合理的，基本已经达到要求。

经过实验和分析可以看出，将本文提出的对NMF算法中施加相关性约束和 稀疏性约束的一种约束NMF算法运用于实际的图像中时，得出的结果较为满意， 具有较强的稀疏性和信号较高的信噪比。