多图正则多核非负矩阵分解高光谱图像解混

1、多图正则多核非负矩阵分解高光谱图像解混高光谱图像1有高的光谱分辨率，包含丰富的图像及光谱信息，而由于高光谱传感器的低空间分辨率和不同纯物质波谱的混合，导致混合像元2的产生，极大地影响了高光谱遥感图像的应用。为改善高光谱图像分解精度，高光谱解混3-5已成为热点，可用线性或非线性方式将混合像元分解，同时提取端元与丰度。端元是混合像元分解出的纯物质光谱，而丰度6则是每个像元中端元的贡献。早期高光谱解混算法主要采用线性混合模型（Linear Mixture Model， LMM），如基于几何的顶点成分分析法（Vertex Component Analysis， VCA）7、最小封闭体积的单纯形法（Mi

2、nimum Volume Enclosing Simplex， MVES）8，基于统计的贝叶斯方法、独立成分分析（Independent Component Analysis， ICA）9和非负矩阵分解（Nonnegative Matrix Factorization， NMF）10算法等。在NMF算法中，先找到一组非负基，然后将原始数据映射到这组基上，且数据在每个基上的表达非负。NMF十分适合应用于高光谱解混，得到端元和丰度矩阵。为保持数据空间固有的流形结构，Yang等人在NMF中加入了图正则算法，称为图正则非负矩阵分解11，该算法利用一种图相似性描述样本之间的关系，充分考虑了局部不变性。由

3、于定义样本之间的关系及其权重矩阵的方法很多，单图的选择至关重要。为解决图选择的问题，Shu等人提出多图正则非负矩阵分解（Multi-graph Regularized Nonnegative Matrix Factorization， MGNMF）12，采用一组已知的多个图拉普拉斯矩阵，通过学习得到的加权参数组合，去逼近原始数据。多图可更准确地刻画样本的相似性，进而更好地表达原始数据的结构。虽然多图解决了图正则非负矩阵分解（Graph Regularized Nonnegative Matrix Factorization， GNMF）算法中图选择困难的问题，但其属于线性解混算法13-14，难

4、以适应真实场景中复杂的非线性光谱混合结构，所以，这一问题也促进了非线性解混算法研究。NMF是一种典型的线性解混算法，不适合提取数据的非线性混合结构15，而核方法可以解决这个问题。核方法16是将非线性关系转变为线性关系的一种过程方法，将低维非线性混合结构的数据映射到高维核空间，在核空间中数据混合结构呈现线性，因此在核空间可实现高光谱图像的非线性解混。Yan等人提出的包含纯像元的核非负矩阵分解（Pure-pixels Kernel Nonnegative Matrix Factorization， pKNMF）与不含纯像元的核非负矩阵分解（Non-pure-pixel Kernel Nonnegt

5、ive Matrix Factorization， npKNMF），是将KNMF17算法直接应用于高光谱数据，得到比NMF好的分类效果。但KNMF的性能很大程度上取决于核函数的选择。多核非负矩阵分解（Multi-kernel Nonnegtive Matrix Factorization， MKNMF）算法采用多个核函数的组合，并为每个核函数设置适当的权重参数。相比KNMF算法，多核NMF可自适应地选择核函数并加权，有更好的映射能力。Yao等人将MKNMF引入图正则NMF算法中，在高光谱数据集中得到了更好的验证18。在许多复杂自然场景中存在大量的非线性混合现象，如沙地和矿物混合区的密集混合现象

6、，以及植被和建筑物覆盖区的多级混合现象。基于线性混合模型的线性解混算法不适合于非线性混合情况，所以，本文提出了一种非线性解混算法多图正则多核非负矩阵分解（Multi-graph Regularized Multi-kernel Nonnegative Matrix Factorization，MGMKNMF），先使用多核函数构造适合于高光谱数据的核空间，然后在多核空间的基础上为目标函数添加多图正则项。本文提出的算法有以下两个优点：（1）与KNMF算法相比，MGMKNMF算法采用核函数权重将多个不同参数的高斯核函数联合起来，并在学习过程中不断更新核函数权重，避免了单核的唯一性，使构造的核空间更合

7、适，也解决了多个高斯核函数权重选择困难的问题。（2）与GNMF和MGNMF算法相比，MGMKNMF算法是非线性方法，在多核空间构造多图，图权重将多个图拉普拉斯矩阵线性组合，并与丰度矩阵最终构成多核空间的多图正则项，且在学习过程中不断更新图权重。相比原空间的单图和多图，多核空间中的多图可更精确地刻画原始数据的非线性流形结构，更适合对真实场景中复杂的非线性光谱混合结构进行非线性解混。本文采用两个真实高光谱数据集Jasper Ridge和Cuprite，并采用广义双线性模型（GBM）19和Hapke20非线性模型分别生成两个模拟数据集，将所提MGMKNMF算法与GNMF、npKNMF、核稀疏非负矩阵

8、分解（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 S

9、quare Error， RMSE）相比其他算法均有较为显著的下降。2 相关工作2.1非负矩阵分解NMFNMF可用于高光谱遥感影像的无监督解混。给定一个数据矩阵X=x1,xNRMN，其中X的每一列都是样本向量。NMF通过矩阵分解将原始高秩矩阵分解为两个低秩矩阵ARMK和SRKN的相乘，并加上非负的约束。XAS（1）非负矩阵分解的目标函数：minf(A,S)=|XAS|2s.t.A0,S0（2）其中：X是原始数据矩阵，在高光谱数据中，A和S分别代表端元矩阵与丰度矩阵。2.2图正则非负矩阵分解GNMFNMF将非负矩阵X分解为基矩阵A和编码矩阵S的乘积(XAS)，用于高光谱图像的无监督解混时，即将高

10、光谱数据集X分解为端元矩阵A与丰度矩阵S的乘积。GNMF将图正则项添加到NMF的目标函数中，改善了未考虑数据流形结构所带来的问题。对有N个样本的数据集X构造一个K近邻图，图中顶点为X中各像素点。Nn为X中样本xn的K近邻集，将每个顶点xn与属于它的Nn连接，并定义图权重矩阵WRNN，原空间中顶点xn和xm间的权重Wnm越大，在子空间中的sn和sm距离也越近。通过权重矩阵W构造图正则项：OG(S,W)=12n,m=1Nsnsm2Wnm=tr(SDST)tr(SWST)=tr(SLST)（3）其中：sn和sm为像元xn和xm在端元基向量上的编码向量，即丰度；Wnm为xn和xm间的权重；S为丰度矩阵

11、；D为对角矩阵，Dnn=Nm=1Wnm；L=DW是图拉普拉斯矩阵；“tr（）”表示矩阵的对角线元素之和。常用定义权重矩阵W的方法有：0-1加权、热核加权和点积加权等。将式（3）与（2）结合，得到GNMF的目标函数：OGNMF(A,S,W)=tr(XXT)2tr(XATST)+tr(ASSTAT)+tr(SLST),（4）其中：为权重参数，GNMF求解约束最小化问题minOGNMF(A,S,W)s.t.A0,S0。2.3多图正则非负矩阵分解MGNMF根据不同定义权重矩阵的模型，可计算相应的图权重矩阵和图拉普拉斯矩阵。MGNMF采用不同数量的最近邻构建图并进行加权，经过算法自动选择，得到最优多图正

12、则项。若已知一组M种模型的图权重矩阵W1，W2，WM和相应的图拉普拉斯矩阵L1，L2，LM，将这M个图权重矩阵线性组合，则相应的图拉普拉斯矩阵也进行相同的线性组合：W=m=1MmWmL=m=1MmLms.t.m=1Mm=1,m0（5）其中，m是第M个图权重矩阵和图拉普拉斯矩阵的权重。MGNMF为一组预先计算得到的候选图确定最佳的图权重，而不是先选最佳图矩阵模型并估计参数。MGNMF的多图正则项为：OMG(S,)=m=1Mmtr(SLmST)s.t.m=1Mm=1,m0（6）多图正则项比单图正则项更精准，GNMF的唯一权重不可靠。MGNMF无需选择唯一的图权重矩阵模型，且通过学习所得图权重向量=

13、1,LT对M个图拉普拉斯矩阵进行最优线性组合。2.4多核非负矩阵分解MKNMFNMF是线性方法，不能很好地处理数据中的非线性结构，而KNMF通过核方法，将原始数据映射到高维核空间：X(X)=(x1),(xn)，可以解决数据的非线性问题。核矩阵K=(X)T(X)，KRNN。则在高维核空间，NMF可以表示为(X)AS。其中A为端元矩阵，S为丰度矩阵。以核空间中样本(X)作为基向量，得到端元矩阵A=(X)F，A=a1,a2,aP，P为端元个数，F矩阵包含核空间中所有样本对构造各端元的贡献，FRNP且ap=n=1Nfnp(xn)，fnp为矩阵F中第n行p列个元素，F矩阵中的第p列为核空间中所有样本对构

14、造第p个端元的贡献。当核函数确定后，KNMF有唯一的核空间，这种选择核空间的方式并不准确。MKNMF算法将L个不同核函数对应的核空间联合起来，以构造一个更合适的希尔伯特空间，这L个不同核空间对应的核矩阵为KlLl=1。核函数权重向量=1,LT将这L个不同的核空间线性组合，组合后的多核核空间的核矩阵为：K=l=1LlKl,s.t.l0,l=1Ll=1（7）通过将学习好的参数带入上式，避免了不同核函数权重分配的问题；而多个核函数可构造出更适合原数据的核空间，比单核更可靠。将式（4）代入KNMF目标函数中，得MKNMF的目标函数：OMKNMF(F,S,)=trl=1LlKl(IFS)(IFS)T（8

15、）3 MGMKNMF本文提出的MGMKNMF在多核空间构造多图，为更新多图，给定参数，用欧几里德距离的平方重新定义多核空间中xn的K近邻集Nn：d(xn,xm)=(xn)(xm)2=K(xn,xn)+K(xm,xm)2K(xn,xm)=l=1LlKl(xn,xn)+Kl(xm,xm)2Kl(xn,xm).（9）多核空间中的多图正则化项为：OMG(S,W)=12n,m=1Nsnsm2Wnm=m=1Mmtr(SLmST),（10）其中，Lm=DmWm是核空间中的图拉普拉斯矩阵。和MGNMF类似，为一组预先计算所得图拉普拉斯矩阵确定最佳权重，精准的多图正则项应用在多核空间将更加合理。最终将式（8）与

16、式（10）结合，得到MGMKNMF的目标函数：minF,S,OMGMKNMF(F,S,)=OMKNMF(F,S,)+OMG(S,W)+2+2=trl=1LlKl(IFS)(IFS)T+m=1Mmtr(SLmST)+2+2=trK(IFS)(IFS)T+tr(SLST)+2+2s.t.F0,S0,l0,l=1Ll=1,m=1Mm=1,m0（11）多核空间可更好地挖掘数据间的非线性关系，在多核空间嵌入多图能更好地表达数据的非线性流形结构。式（11）中的2可防止参数过度偏向到一个核函数中；2项可避免偏向到一个权重构造函数中；，均是权衡上式所用的权衡参数，其值均为非负。因是约束多图正则项的参数，与约束

17、核权重和图权重的，相比，应大于，。根据参考文献21和多次实验的结果，本文实验中：Cuprite数据=100，Jasper Ridge数据=20，HAPKE模拟数据=20，GBM模拟数据=20；所有数据的均为10，均为10。根据参考文献21和22，式（11）中的MGMKNMF目标函数是非凸的，无法得到全局最小值。MGMKNMF采用分步迭代策略优化目标函数，可得到局部最小值21-22，具体如下：（1）固定和，更新F，S。式（11）可写为：minF,StrK(IFS)(IFS)T+tr(SLST)s.t.F0,S0.(12)假设lp和pn分别是flp和spn的拉格朗日乘子，令lp=，pn=，则式（1

18、2）的拉格朗日函数为：=minF,StrK(IFS)(IFS)T+tr(SLST)+tr(FT)+tr(ST).（13）式（13）分别对F和S求偏导得：F=KFSSTKST+S=STFTKFKF+S(L+(L)T)+（14）结合卡罗需-库恩-塔克（KKT）条件：lpflp=0，pnspn=0有：(KFSST)flp(KST)flp+flp=0(STFTKF)spn(KF)spn+(S(L+(L)T)spn+spn=0（15）得到F和S的乘法更新法则分别为：flp(KST)lp(KFSST)flp（16）spn(FTK+SW)(FTKFS+SD)spn（17）（2）固定F，S和，更新，式（11）

19、可写为：mintrl=1LlKl(IFS)(IFS)T+2=trl=1LlKlZZT+2=l=1Llgl+l=1Ll2s.t.l0,l=1Ll=1（18）其中，Z=IFS，gl=trKlZZT。（3）固定F和S和，更新，式（11）可以写为：minm=1Mmtr(SLmST)+2s.t.m=1Mm=1,m0（19）式（18）和式（19）的约束二次规划问题可根据文献23中的方法来解决。MGMKNMF解混算法总结如下：MGMKNMF解混算法输入：原始高光谱数据X，L个核矩阵K1,Kl,KLK1,Kl,KL，M个图拉普拉斯矩阵L1,Lm,LML1,Lm,LM，最大迭代次数TStep1.初始化矩阵F0和

20、S0，Step2.初始化核权重变量0l=1/L,l=1,Ll0=1/L,l=1,L，初始化图权重0m=1/M,m=1,Mm0=1/M,m=1,Mfort=1 toT do通过式（7）更新图GtGt和相应的拉普拉斯矩阵LtLt；通过式（11）和式（12）更新FtFt和StSt；通过式（13）和式（14）更新核权重tt和图权重tt；end输出：F=Ft1,S=St1,F=Ft-1,S=St-1,=t1,=t1=t-1,=t-14 实验结果与分析将所提MGMKNMF算法与GNMF、npKNMF、KSNMF、KDP和MGKNMF解混算法进行对比，使用SAD和RMSE作为评估指标，采用HAPKE和GBM

21、模拟数据，以及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|sisi|2)1/2，其中

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

23、mly generated by spectrum library4.3模拟数据实验结果及分析表14分别显示了端元数目为6，SNR不同时，各算法在HAPKE和GBM模拟数据上的SAD、RMSE值。从表中可以看出，随着SNR值的增加，SAD和RMSE的值均呈下降趋势，而 MGMKNMF算法与其他6种算法相比具有较好的准确度。NMF和GNMF因为属于线性解混算法，对非线性数据解混精度都很差，但GNMF算法在NMF的基础上增加了图正则约束，揭示了数据内在固有的流形结构，相比NMF算法有进一步提升；剩下的5种算法均利用了核函数的概念，更适合于非线性数据。npKNMF算法在核空间中应用了NMF，相比NM

24、F准确度有所提高；而KSNMF在丰度中加了L1范数使丰度更加稀疏，得到了更好的准确率；KDP在进行端元选择时使用了大型光谱库，相比其余两者更加准确，MGKNMF在单核空间中构造多图正则项。而提出的MGMKNMF算法不仅用多图来刻画数据内在流形结构，更使用多核学习找到合适的核参数与核函数，进一步增加了算法的准确度。表1不同SNR下HAPKE模型各算法的SAD值Tab.1SAD value of each algorithm of HAPKE model under different SNRSNR/dBNMFGNMFnpKNMFKSNMFKDPMGKNMFMGMKNMF100.494 80.45

25、6 70.170 20.169 30.167 10.190 10.213 3200.382 90.309 40.164 80.160 40.154 30.157 60.154 1300.341 70.296 00.162 10.153 10.148 20.129 4.0.124 8400.308 40.251 90.159 60.149 20.142 70.146 60.137 5表2不同SNR下GBM模型各算法的SAD值Tab.2SAD value of each algorithm of GBM model under different SNRSNR/dBNMFGNMFnpKNMFKSN

26、MFKDPMGKNMFMGMKNMF100.401 10.283 90.180 30.172 30.171 40.169 30.258 8200.395 90.236 40.174 20.166 80.164 30.145 50.121 9300.341 30.214 20.168 10.161 10.160 80.139 00.113 1400.300 10.199 80.161 30.159 10.157 70.149 90.145 6表3不同SNR下HAPKE模型各算法的RMSE值Tab.3RMSE value of each algorithm of HAPKE model under

27、 different SNRSNR/dBNMFGNMFnpKNMFKSNMFKDPMGKNMFMGMKNMF100.185 00.186 90.079 20.068 90.067 20.099 00.091 0200.162 40.183 90.076 90.063 70.061 10.061 80.060 1300.149 50.180 10.071 30.058 70.058 00.052 30.049 0400.129 90.178 30.066 20.052 30.051 80.051 50.050 6表4不同SNR下GBM模型各算法的RMSE值Tab.4RMSE value of e

28、ach algorithm of GBM model under different SNRSNR/dBNMFGNMFnpKNMFKSNMFKDPMKGNMFMGMKNMF100.321 70.290 60.081 70.071 10.070 80.079 00.112 3200.307 40.253 10.079 40.062 40.061 80.060 60.059 7300.292 30.233 70.073 60.060 30.059 30.052 20.050 1400.258 60.201 20.069 10.058 40.058 00.057 70.057 0表58分别显示了SN

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

30、orithm of HAPKE model under different number of endmembersEndmemberNMFGNMFnpKNMFKSNMFKDPMGKNMFMGMKNMFP=60.308 40.281 70.159 60.149 00.142 20.145 70.141 1P=50.251 00.218 70.124 00.118 20.107 70.096 80.081 4P=40.198 50.195 40.100 40.089 80.081 20.080 20.078 9P=30.172 90.143 80.073 50.068 20.060 00.044

31、 30.034 3表6不同端元数目下GBM模型各算法的SAD值Tab.6SAD value of each algorithm of GBM model under different number of endmembersEndmemberNMFGNMFnpKNMFKSNMFKDPMGKNMFMGMKNMFP=60.271 10.278 90.169 50.159 80.158 00.148 30.147 1P=50.268 30.241 90.125 10.121 50.113 60.100 20.095 1P=40.209 20.202 10.091 40.085 40.081 70.

32、070 00.063 2P=30.182 00.143 50.067 90.060 50.057 20.042 30.029 9表7不同端元数目下HAPKE模型各算法的RMSE值Tab.7RMSE value of each algorithm of HAPKE model under different number of endmembersEndmemberNMFGNMFnpKNMFKSNMFKDPMGKNMFMGMKNMFP=60.189 90.176 40.084 10.068 10.066 90.060 30.059 4P=50.175 30.162 10.078 30.062 3

33、0.061 10.060 10.056 3P=40.151 50.160 10.074 50.058 90.057 40.054 70.052 3P=30.203 00.157 70.066 20.055 30.051 80.051 40.050 8表8不同端元数目下GBM模型各算法的RMSE值Tab.8RMSE value of each algorithm of GBM model under different number of endmembersEndmemberNMFGNMFnpKNMFKSNMFKDPMGKNMFMGMKNMFP=60.227 20.200 10.086 30.

34、079 20.077 50.060 10.061 0P=50.180 40.173 40.081 20.068 30.066 20，0 6220.060 3P=40.153 40.149 70.079 90.065 10.064 30.064 00.063 3P=30.119 60.139 60.069 10.058 40.058 00.041 20.035 44.4真实数据实验结果及分析本文采用真实地物Cuprite和Jasper Ridge数据集对MGMKNMF的有效性进行验证，两个数据集均可以在https：/rslab.ut.ac.ir/data中下载。Cuprite是高光谱解混研究常用

35、的数据集，包含美国内华州Cuprite矿区。在除去低信噪比和吸水通道后，留有188个通道可以使用。每张图像大小为250191，共有12种类别。表9总结了在Cuprite数据集各类算法的SAD值。可以看出，MGMKNMF算法的平均SAD值是最优的。表9不同算法在Cuprite数据的SAD值Tab.9SAD values of Cuprite data by different algorithmsItemNMFGNMFnpKNMFKSNMFKDPMGKNMFMGMKNMFAlunite0.297 90.243 40.063 20.085 90.062 90.078 70.069 5Andradi

36、te0.383 90.262 60.079 50.101 40.204 50.066 00.079 4Buddingtonite0.367 20.266 60.079 60.073 70.116 10.088 50.119 6Dumortierite0.283 90.574 30.156 20.071 70.071 20.084 80.082 9Kaolinite_10.235 90.336 20.080 20.061 30.081 30.082 10.087 5Kaolinite_20.309 50.344 30.262 20.189 00.114 00.056 50.073 3Muscov

37、ite0.379 40.318 30.113 70.052 10.147 20.102 40.104 1Montmorillonite0.473 10.355 70.131 00.140 20.107 00.054 50.055 6Nontronite0.451 60.451 60.129 30.064 70.078 80.121 40.104 3Pyrope0.292 70.455 00.054 10.142 00.074 50.117 80.065 5Sphene0.282 10.287 90.091 30.279 80.097 60.201 50.197 1Chalcedony0.318 60.318 60.213 20.149 50.132 70.071 50.065 9Average0.339 60.337 10.120 80.117 60.107 80.093 80.092 1图2为MGMKNMF算法在Cuprite数据集上的丰度图。结合表9，在此算法下，12种端元对应的丰度图应是最优的，可清晰地看出各类别的划分情况。图2MGM