刘 敬，李康欣，张 悠，刘 逸

（1. 西安邮电大学 电子工程学院，陕西 西安 710121；2. 西安电子科技大学 电子工程学院，陕西 西安 710071）

1 引言

高光谱图像［1］有高的光谱分辨率，包含丰富的图像及光谱信息，而由于高光谱传感器的低空间分辨率和不同纯物质波谱的混合，导致混合像元［2］的产生，极大地影响了高光谱遥感图像的应用。为改善高光谱图像分解精度，高光谱解混［3-5］已成为热点，可用线性或非线性方式将混合像元分解，同时提取端元与丰度。端元是混合像元分解出的纯物质光谱，而丰度［6］则是每个像元中端元的贡献。

早期高光谱解混算法主要采用线性混合模型（Linear Mixture Model，LMM），如基于几何的顶点成分分析法（Vertex Component Analysis，VCA）［7］、最小封闭体积的单纯形法（Minimum Volume Enclosing Simplex，MVES）［8］，基于统计的贝叶斯方法、独立成分分析（Independent Component Analysis，ICA）［9］和非负矩阵分解（Nonnegative Matrix Factorization，NMF）［10］算法等。在NMF 算法中，先找到一组非负基，然后将原始数据映射到这组基上，且数据在每个基上的表达非负。NMF 十分适合应用于高光谱解混，得到端元和丰度矩阵。为保持数据空间固有的流形结构，Yang 等人在NMF 中加入了图正则算法，称为图正则非负矩阵分解［11］，该算法利用一种图相似性描述样本之间的关系，充分考虑了局部不变性。由于定义样本之间的关系及其权重矩阵的方法很多，单图的选择至关重要。为解决图选择的问题，Shu 等人提出多图正则非负矩阵分解（Multi-graph Regularized Nonnegative Matrix Factorization，MGNMF）［12］，采用一组已知的多个图拉普拉斯矩阵，通过学习得到的加权参数组合，去逼近原始数据。多图可更准确地刻画样本的相似性，进而更好地表达原始数据的结构。虽然多图解决了图正则非负矩阵分解（Graph Regularized Nonnegative Matrix Factorization，GNMF）算法中图选择困难的问题，但其属于线性解混算法［13-14］，难以适应真实场景中复杂的非线性光谱混合结构，所以，这一问题也促进了非线性解混算法研究。

NMF 是一种典型的线性解混算法，不适合提取数据的非线性混合结构［15］，而核方法可以解决这个问题。核方法［16］是将非线性关系转变为线性关系的一种过程方法，将低维非线性混合结构的数据映射到高维核空间，在核空间中数据混合结构呈现线性，因此在核空间可实现高光谱图像的非线性解混。Yan 等人提出的包含纯像元的核非负矩阵分解（Pure-pixels Kernel Nonnegative Matrix Factorization，pKNMF）与不含纯像元的核非负矩阵分解（Non-pure-pixel Kernel Nonnegtive Matrix Factorization，npKNMF），是将KNMF［17］算法直接应用于高光谱数据，得到比NMF 好的分类效果。但KNMF 的性能很大程度上取决于核函数的选择。多核非负矩阵分解（Multi-kernel Nonnegtive Matrix Factorization，MKNMF）算法采用多个核函数的组合，并为每个核函数设置适当的权重参数。相比KNMF 算法，多核NMF 可自适应地选择核函数并加权，有更好的映射能力。Yao 等人将MKNMF 引入图正则NMF 算法中，在高光谱数据集中得到了更好的验证［18］。

在许多复杂自然场景中存在大量的非线性混合现象，如沙地和矿物混合区的密集混合现象，以及植被和建筑物覆盖区的多级混合现象。基于线性混合模型的线性解混算法不适合于非线性混合情况，所以，本文提出了一种非线性解混算法——多图正则多核非负矩阵分解（Multigraph Regularized Multi-kernel Nonnegative Matrix Factorization，MGMKNMF），先使用多核函数构造适合于高光谱数据的核空间，然后在多核空间的基础上为目标函数添加多图正则项。本文提出的算法有以下两个优点：（1）与KNMF 算法相比，MGMKNMF 算法采用核函数权重将多个不同参数的高斯核函数联合起来，并在学习过程中不断更新核函数权重，避免了单核的唯一性，使构造的核空间更合适，也解决了多个高斯核函数权重选择困难的问题。（2）与GNMF 和MGNMF 算 法 相 比，MGMKNMF 算 法 是 非 线 性方法，在多核空间构造多图，图权重将多个图拉普拉斯矩阵线性组合，并与丰度矩阵最终构成多核空间的多图正则项，且在学习过程中不断更新图权重。相比原空间的单图和多图，多核空间中的多图可更精确地刻画原始数据的非线性流形结构，更适合对真实场景中复杂的非线性光谱混合结构进行非线性解混。本文采用两个真实高光谱数据集Jasper Ridge 和Cuprite，并采用广义双线性模型（GBM）［19］和Hapke［20］非线性模型分别生成两个模拟数据集，将所提MGMKNMF 算法与GNMF、npKNMF、核稀疏非负矩阵分解（Kernel Sparse Nonnegative Matrix Factorization，KSNMF）、基于核的字典剪枝非线性光谱解混（Kernel-based Nonlinear Spectral Unmixing with Dictionary Pruning，KDP）、多图正则核非负矩阵（Multi-graph Regularized Kernel Nonnegative Matrix Factorization， MGKNMF）算法比较，实验结果表明，MGMKNMF的光谱角距离（Spectral Angel Distance，SAD）和 均 方 根 误 差（Root Mean Square Error，RMSE）相比其他算法均有较为显著的下降。

2 相关工作

2.1 非负矩阵分解NMF

NMF 可用于高光谱遥感影像的无监督解混。给定一个数据矩阵X=[x1，…，xN]∈RM×N，其中X的每一列都是样本向量。NMF 通过矩阵分解将原始高秩矩阵分解为两个低秩矩阵A∈RM×K和S∈RK×N的相乘，并加上非负的约束。

非负矩阵分解的目标函数：

其中：X是原始数据矩阵，在高光谱数据中，A和S分别代表端元矩阵与丰度矩阵。

2.2 图正则非负矩阵分解GNMF

NMF 将非负矩阵X分解为基矩阵A和编码矩阵S的乘积(X≈AS)，用于高光谱图像的无监督解混时，即将高光谱数据集X分解为端元矩阵A与丰度矩阵S的乘积。GNMF 将图正则项添加到NMF 的目标函数中，改善了未考虑数据流形结构所带来的问题。对有N个样本的数据集X构造一个K近邻图，图中顶点为X中各像素点。Nn为X中样本xn的K近邻集，将每个顶点xn与属于它的Nn连接，并定义图权重矩阵W∈RN×N，原空间中顶点xn和xm间的权重Wnm越大，在子空间中的sn和sm距离也越近。通过权重矩阵W构造图正则项：

其中：sn和sm为像元xn和xm在端元基向量上的编码向量，即丰度；Wnm为xn和xm间的权重；S为丰度矩阵；D为对角矩阵，D-W是图拉普拉斯矩阵；“tr（·）”表示矩阵的对角线元素之和。常用定义权重矩阵W的方法有：0-1 加权、热核加权和点积加权等。

将式（3）与（2）结合，得到GNMF 的目标函数：

其中：α为权重参数，GNMF 求解约束最小化问题minOGNMF(A，S，W) s.t.A≥0，S≥0。

2.3 多图正则非负矩阵分解MGNMF

根据不同定义权重矩阵的模型，可计算相应的图权重矩阵和图拉普拉斯矩阵。MGNMF 采用不同数量的最近邻构建图并进行加权，经过算法自动选择，得到最优多图正则项。若已知一组M种模型的图权重矩阵｛W1，W2，…，WM｝和相应的图拉普拉斯矩阵｛L1，L2，…，LM｝，将这M个图权重矩阵线性组合，则相应的图拉普拉斯矩阵也进行相同的线性组合：

其中，γm是第M个图权重矩阵和图拉普拉斯矩阵的权重。MGNMF 为一组预先计算得到的候选图确定最佳的图权重，而不是先选最佳图矩阵模型并估计参数。MGNMF 的多图正则项为：

多图正则项比单图正则项更精准，GNMF 的唯一权重不可靠。MGNMF 无需选择唯一的图权重矩阵模型，且通过学习所得图权重向量γ=[γ1，…，γL]T对M个图拉普拉斯矩阵进行最优线性组合。

2.4 多核非负矩阵分解MKNMF

NMF 是线性方法，不能很好地处理数据中的非线性结构，而KNMF 通过核方法，将原始数据 映 射 到 高 维 核 空 间 ：X→Φ(X)=[Φ(x1)，…，Φ(xn)]，可以解决数据的非线性问题。核矩阵K=Φ(X)TΦ(X)，K∈RN×N。则在高维核空间，NMF 可以表示为Φ(X)≈AS。其中A为端元矩阵，S为丰度矩阵。以核空间中样本Φ(X) 作为基向量，得到端元矩阵A=Φ(X)F，A=[a1，a2，…，aP]，P为端元个数，F矩阵包含核空间中所有样本对构造各端元的贡献，F∈RN×P且为矩阵F中 第n行p列个元素，F矩阵中的第p列为核空间中所有样本对构造第p个端元的贡献。当核函数确定后，KNMF 有唯一的核空间，这种选择核空间的方式并不准确。

MKNMF 算法将L个不同核函数对应的核空间联合起来，以构造一个更合适的希尔伯特空间，这L个不同核空间对应的核矩阵为核函数权重向量τ=[τ1，…，τL]T将这L个不同的核空间线性组合，组合后的多核核空间的核矩阵为：

通过将学习好的参数τ带入上式，避免了不同核函数权重分配的问题；而多个核函数可构造出更适合原数据的核空间，比单核更可靠。将式（4）代入KNMF 目标函数中，得MKNMF 的目标函数：

3 MGMKNMF

本文提出的MGMKNMF 在多核空间构造多图，为更新多图，给定参数τ，用欧几里德距离的平方重新定义多核空间中xn的K近邻集Nτ n：

多核空间中的多图正则化项为：

最终将式（8）与式（10）结合，得到MGMKNMF 的目标函数：

多核空间可更好地挖掘数据间的非线性关系，在多核空间嵌入多图能更好地表达数据的非线性流形结构。式（11）中的‖τ‖2可防止参数τ过度偏向到一个核函数中；‖γ‖2项可避免γ偏向到一个权重构造函数中；α，β，μ均是权衡上式所用的权衡参数，其值均为非负。因α是约束多图正则项的参数，与约束核权重和图权重的β，μ相比，α应大于β，μ。根据参考文献［21］和多次实验的结果，本文实验中：Cuprite 数据α=100，Jasper Ridge 数据α=20，HAPKE 模拟数据α=20，GBM 模拟数据α=20；所有数据的β均为10，μ均为10。

根据参考文献［21］和［22］，式（11）中的MGMKNMF 目标函数是非凸的，无法得到全局最小值。MGMKNMF 采用分步迭代策略优化目标函数，可得到局部最小值［21-22］，具体如下：

（1）固定τ和γ，更新F，S。式（11）可写为：

假设φlp和ψpn分别是flp和spn的拉格朗日乘子，令φlp=φ，ψpn=ψ，则式（12）的拉格朗日函数为：

式（13）分别对F和S求偏导得：

结合卡罗需-库恩-塔克（KKT）条件：φlp flp=0，ψpn spn=0 有：

得到F和S的乘法更新法则分别为：

（2）固定F，S和γ，更新τ，式（11）可写为：

其中，Z=I-FS，gl=tr[Kl ZZT]。

（3）固 定F和S和τ，更 新γ，式（11）可 以写为：

式（18）和式（19）的约束二次规划问题可根据文献［23］中的方法来解决。

MGMKNMF 解混算法总结如下：

MGMKNMF 解混算法输入：原始高光谱数据X，L 个核矩阵{ K1，…，Kl，…，KL}，M 个图拉普拉斯矩阵{ L1，…，Lm，…，LM}，最大迭代次数T Step1. 初始化矩阵F0和S0，Step2. 初始化核权重变量τ0 l =1 L，l=1，…，L，初始化图权重γ0m=1 M，m=1，…，M for t=1 to T do通过式（7）更新图Gτt 和相应的拉普拉斯矩阵Lτt；通过式（11）和式（12）更新Ft 和St；通过式（13）和式（14）更新核权重τt 和图权重γt；end输出：F=Ft-1，S=St-1，τ=τt-1，γ=γt-1

4 实验结果与分析

将 所 提 MGMKNMF 算 法 与 GNMF、npKNMF、KSNMF、KDP 和MGKNMF 解 混 算法进行对比，使用SAD 和RMSE 作为评估指标，采用HAPKE 和GBM 模拟数 据，以及Cuprite 和Jasper Ridge 真实数据验证该算法的有效性。所有实验中，我们选择0-1 加权图、热核加权图和点积加权图的图权重矩阵模型构成最终图正则项。多核函数选择高斯核函数，所选取的核参数为1/32，1/16，1/8，1/4，1/2，1，2，4，8，16，32 以构成不同的核函数。各算法的迭代次数T均设置为200。

4.1 评价标准

SAD值反映了解混所得端元光谱与原端元光谱之间相似性，定义为VSAD=arccos(aTb/||a||||b||)，其中，a，b是两个端元光谱。

RMSE 值反映解混所得丰度与实验室测量的实际丰度间的差别，定义为VRMSEsi=(1/N|si-s^i|2)1/2，其中，si和s^i分别是实际丰度和解混所得丰度。

4.2 模拟数据实验与分析

本文采用HAPKE 和GBM 两种非线性模型生成模拟数据集。HAPKE 模型是一种紧密混合模型；GBM 是一种双线性混合模型，是LMM 线性结构与端元间的二次散射项的加权组合。本文从美国光谱库（USGS）随机选择6 种地物光谱作为端元，如图1。并且丰度矢量满足丰度非负和丰度和为一的约束。最后模拟生成的高光谱每个大小为20×20，每个像素波段为224，并加入不同信噪比（Signal to Noise Ratio，SNR）的零均值高斯噪声来更好地接近真实数据。

图1 光谱库中随机生成的端元光谱Fig.1 Endmember spectra randomly generated by spectrum library

4.3 模拟数据实验结果及分析

表1～4 分别显示了端元数目为6，SNR 不同时，各算法在HAPKE 和GBM 模拟数据上的SAD、RMSE 值。从表中可以看出，随着SNR 值的增加，SAD 和RMSE 的值均呈下降趋势，而MGMKNMF 算法与其他6 种算法相比具有较好的准确度。NMF 和GNMF 因为属于线性解混算法，对非线性数据解混精度都很差，但GNMF算法在NMF 的基础上增加了图正则约束，揭示了数据内在固有的流形结构，相比NMF 算法有进一步提升；剩下的5 种算法均利用了核函数的概念，更适合于非线性数据。npKNMF 算法在核空间中应用了NMF，相比NMF 准确度有所提高；而KSNMF 在丰度中加了L1范数使丰度更加稀疏，得到了更好的准确率；KDP 在进行端元选择时使用了大型光谱库，相比其余两者更加准确，MGKNMF 在单核空间中构造多图正则项。而提出的MGMKNMF 算法不仅用多图来刻画数据内在流形结构，更使用多核学习找到合适的核参数与核函数，进一步增加了算法的准确度。

表5～8 分别显示了SNR 为40 dB，端元数目P不同时，各算法在这2 个模型的模拟数据上的SAD 与RMSE 值。因每次生成模拟数据时是在USGS 库随机选取端元，导致表5～8 的结果在P=6 时，与表1～4 中SNR 为40 dB 时的结果不同。可以看出，因为模拟数据非线性程度高，端元数的增加会导致算法解混能力的下降。相比NMF 算 法，MGMKNMF 算 法 在HAPKE 数 据中，端元为6 时，SAD 值和RMSE 值分别减少了约0.17 和0.13；在GBM 数据中，分别减少了约0.12 和0.17。和其余算法相比，MGMKNMF 算法也基本保持着最优结果。

表1 不同SNR 下HAPKE 模型各算法的SAD 值Tab.1 SAD value of each algorithm of HAPKE model under different SNR

表2 不同SNR 下GBM 模型各算法的SAD 值Tab.2 SAD value of each algorithm of GBM model under different SNR

表3 不同SNR 下HAPKE 模型各算法的RMSE 值Tab.3 RMSE value of each algorithm of HAPKE model under different SNR

表4 不同SNR 下GBM 模型各算法的RMSE 值Tab.4 RMSE value of each algorithm of GBM model under different SNR

表5 不同端元数目下HAPKE 模型各算法的SAD 值Tab.5 SAD value of each algorithm of HAPKE model under different number of endmembers

表6 不同端元数目下GBM 模型各算法的SAD 值Tab.6 SAD value of each algorithm of GBM model under different number of endmembers

表7 不同端元数目下HAPKE 模型各算法的RMSE 值Tab.7 RMSE value of each algorithm of HAPKE model under different number of endmembers

表8 不同端元数目下GBM 模型各算法的RMSE 值Tab.8 RMSE value of each algorithm of GBM model under different number of endmembers

4.4 真实数据实验结果及分析

本文采用真实地物Cuprite 和Jasper Ridge 数据集对MGMKNMF 的有效性进行验证，两个数据 集 均 可 以 在https：//rslab. ut. ac. ir/data 中下载。

Cuprite 是高光谱解混研究常用的数据集，包含美国内华州Cuprite 矿区。在除去低信噪比和吸水通道后，留有188 个通道可以使用。每张图像大小为250×191，共有12 种类别。

表9 总结了在Cuprite 数据集各类算法的SAD 值。可以看出，MGMKNMF 算法的平均SAD 值是最优的。

图2 为MGMKNMF 算法在Cuprite 数据集上的丰度图。结合表9，在此算法下，12 种端元对应的丰度图应是最优的，可清晰地看出各类别的划分情况。

图2 MGMKNMF 算法在Cuprite 数据上的丰度图Fig.2 Abundance graphs of MGMKNMF algorithm on Cuprite

表9 不同算法在Cuprite 数据的SAD 值Tab.9 SAD values of Cuprite data by different algorithms

由于Jasper Ridge 数据集太复杂无法处理，我们仅考虑100×100 像素的子图像，子图像的第一个像素从原始图像中的第（105，269）像素开始，且去除了低噪声和水吸收波段，保留198 个通道。所以，该高光谱遥感数据中只有树、水、土壤和道路这4 类端元。

从表10 可以看到，提出的MGMKNMF 算法在Jasper Ridge 数据集上仍然有效，平均SAD 值依旧是最优的。KDP 算法在KSNMF算法的基础上利用光谱库挑选端元，而MGKNMF 在单核空间应用多图，精度仅次于MGMKNMF。

表10 不同算法在Jasper Ridge 数据的SAD 值Tab.10 SAD values of Jasper Ridge data by different algorithms

图3 为Jasper Ridge 数据集各算法解混出的丰度图，从左至右分别对应树、水、土壤和道路这4 个端元。

图3 各算法在Jasper Ridge 的丰度图Fig.3 Abundance graphs of each algorithm on Jasper Ridge

从图3 可看出，NMF 算法与真实地物丰度相差很大，尤其在水这一端元上；GNMF 算法比NMF 算法的丰度图清晰。而其余的5种算法均在核空间进行，丰度图明显优于NMF 和GNMF，其中水的丰度图更接近真实地物。相比其他算法，在核空间进行的算法尽管其丰度图的区别甚微，但结合表10 的数据，可 以 看 到，MKNGNMF 算 法 的SAD 值 是最优的。

5 结论

本文提出了一种MGMKNMF 高光谱非线性解混算法。该算法用多个核函数构造出了多核空间，且在学习过程中不断更新核函数权重，更有利于揭示原始数据的非线性结构；并在多核空间用图权重向量将多个图拉普拉斯矩阵线性组合，与丰度矩阵最终构成多核空间的多图正则项。相比单图正则项，多核空间的多图正则项将更加逼近原始数据的非线性流形结构。基于2 个真实数据集和2 个模拟数据集的实验结果表明，相比GNMF、npKNMF、KSNMF、MGKNMF 和KDP 算法，MGMKNMF 算法确实是最优的，它更合适于复杂的高光谱数据。