一种基于最小距离和稀疏图正则约束的非负矩阵解混算法

A Nonnegative Matrix Unmixing Algorithm Based on Minimum Distance and Sparse Graph Regularity

摘要: 针对非负矩阵解混模型存在容易陷入局部极小值、受噪声影响较大以及混合像元中存在“异物同谱”端元的情况下解混精度较低等问题，为此，本文提出了一种基于最小距离和稀疏图正则约束的非负矩阵分解算法(DLGNMF)。该算法将空间信息与光谱信息相结合，在目标函数中引入单形体距离最小化约束作为端元的约束条件，弥补解混过程中高光谱数据空间几何特性的缺失，将稀疏约束和图正则化约束作为丰度的约束条件，保证了高光谱图像中端元分布的全局稀疏性与局部相似性，更好的改善了非负矩阵模型存在的问题，提高了解混的精度。

Abstract: Aiming at the problems of non-negative matrix unmixing model, such as easy to fall into the local minimum, being greatly affected by noise, and low unmixing accuracy in the case of “foreign matter having the same spectrum” end elements in mixed pixels, etc., in this paper, a non-negative matrix factorization algorithm (DLGNMF) based on minimum distance and regular constraints of sparse graphs is proposed. The algorithm combines spatial information with spectral information. In the objective function, the constraint of minimizing the distance of a single body is introduced as the constraint condition of the end-member. The sparsity constraint and graph regularization constraint are taken as the constraints of abundance, which ensure the global sparsity and local similarity of the end-member distribution in the hyperspectral image, solve the problems of non-negative matrix model and improve the accuracy of solution unmixing.

1. 引言

高光谱图像其光谱分辨率较高，在资源调查、矿床勘探、灾情预警、和军事伪装识别等领域具有重要的应用价值 [1]。由于传感器平台的局限性以及地物真实分布的复杂性，存在着多种地物组合而成的混合像元，限制了高光谱图像应用的发展 [2]。为此，需从高光谱图像中提取出组成混合像元的纯物质光谱(端元)及其所占的比例(丰度)，这个过程称之为解混 [3]。

Figure 1. Algorithm flow chart

图1. 算法流程图

在解混的过程中，端元的提取决定丰度反演的结果，最终会影响光谱解混的准确性。针对端元提取，基于纯像元模型、基于最小体积模型假设高光谱图像中存在纯像元 [4]，这与实际情况不一致，因此将目光转向不需要假定纯像元存在的基于统计模型的方法 [5]。Lee等人提出一种经典的基于非负矩阵的盲源分解算法，即将任意一个非负矩阵分解为两个低秩非负矩阵的乘积 [6]。

但是，由于非负矩阵分解算法存在非凸性，导致容易陷入局部极小值，以及当混合像元中存在“异物同谱”端元的情况下解混精度较低的问题，为此，本文在非负矩阵模型中引入单形体距离最小化约束作为端元的约束条件，将稀疏约束和图正则化约束作为丰度的约束条件，提出了一种基于最小距离和稀疏图正则的非负矩阵解混算法(DLGNMF)，旨在解决上述问题。流程图如图1所示。

2. 理论介绍

2.1. 线性混合的高光谱非负矩阵解混模型

目前应用较多的混合像元分解模型为线性光谱混合模型 [4] [5]。该模型假设高光谱图像由各个端元线性混合而成。

$X=AS+n$ (1)

其中，丰度服从两个约束条件：1) 非负性约束；2) 和为1约束。X^L×N为高光谱图像；A^L×P为端元矩阵；S^P×N为丰度矩阵；n为误差影响；L为波段数；P为端元数量；N为图像中所含像元总数。

该模型求解的目标函数如下：

$\min f\left(A,S\right)=\frac{1}{2}{\Vert X-AS\Vert }_F^2$ (2)

采用乘性迭代更新法则来对NMF进行求解，得到最终更新法则为：

$A\leftarrow A\frac{X{S}^{\text{T}}}{AS{S}^{\text{T}}}$ (3)

$S\leftarrow S\frac{A^{\text{T}}X}{A^{\text{T}}AS}$ (4)

2.2. 基于稀疏图正则的非负矩阵解混算法

真实的高光谱图像数据中，大多数的像元仅仅包含少量的端元成分，并不是每个像元都包含全部端元成分，即丰度矩阵具有一定的稀疏性 [7]。如果两个像元X_i、X_j离的很近，则两个像元具有一定的相似性，其丰度值S_i与S_j也具有一定的相似性，该假设称之为流行正则假设 [8]，即丰度矩阵具有局部相似性。将稀疏性与流行正则假设作为约束条件，添加到非负矩阵解混模型中 [9]。其目标函数为：

$f\left(A,S\right)=\frac{1}{2}{\Vert X-AS\Vert }_F^2+\lambda {\Vert S\Vert }_{1/2}+\frac{\mu }{2}Tr\left(SL{S}^{\text{T}}\right)$ (5)

其中 $\lambda$ 为稀疏系数； $\mu$ 为图正则化参数； ${\Vert \cdot \Vert }_{1/2}$ 表示矩阵的1/2范数；L为拉普拉斯矩阵； $Tr(\cdot )$ 表示矩阵的迹。

但是GLNMF算法忽略了高光谱图像空间几何性的特点，针对混合像元中存在“异物同谱”端元的情况下，解混精度并没有提升较多，并且对噪声较敏感，在信噪比较低的情况下，解混结果精度偏低。

3. 基于最小距离和稀疏图正则的非负矩阵解混算法

3.1. 算法描述

根据高光谱图像的空间几何性可知，由端元所构成的单形体包含图像中全部像元点，其单形体体积越小，像元点在其中的分布越合理。由于单形体体积计算过于复杂，因此将单形体各个顶点到单形体中心点的距离和最小来代替单形体体积最小。这样既可以降低算法的复杂度，提高解混的效率，又可以使全部像元点在空间中分布更加合理，更好的提升混合像元中“异物同谱”端元存在情况下的解混精度。DLGNMF算法目标函数为：

$f\left(A,S\right)=\frac{1}{2}{\Vert X-AS\Vert }_F^2+\frac{\alpha }{2}Tr\left(AQ{A}^{\text{T}}\right)+\lambda {\Vert S\Vert }_{1/2}+\frac{\mu }{2}Tr\left(SL{S}^{\text{T}}\right)$ (6)

$Tr\left(AQ{A}^{\text{T}}\right)$ 表示最小距离约束， $\alpha$ 表示控制端元约束条件的平衡参数。

将ED定义为单形体各个顶点(即端元)到所有端元中心点的距离和，ED [2] 和Q [2] 的表达式如式(7)、(8)和(9)所示：

$ED\left(A\right)={\displaystyle \underset{i=1}{\overset{P}{\sum }}{\Vert a_i-\bar{a}\Vert }_2^2},\bar{a}=\frac{1}{P}{\displaystyle \underset{i=1}{\overset{P}{\sum }}{a}_i}$ (7)

$ED\left(A\right)=Tr\left(Q^{\text{T}}{A}^{\text{T}}AQ\right)=Tr\left(AQ{A}^{\text{T}}\right)$ (8)

$Q=I-\frac{1}{P}{I}_P{I}_P^{\text{T}}$ (9)

其中 $a_i\left(i=1,2,\cdots ,P\right)$ 为第i个端元；I_P是P × P的单位矩阵；

${\Vert S\Vert }_{1/2}={\displaystyle \underset{i=1}{\overset{P}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{S}_{ij}}}$ (10)

为了描述数据的内在几何结构，采用热核权重函数来构造权重矩阵W，进而构造像元的最近邻图来描述流行正则假设。其定义如下：

$W_{i,j}=\{\begin{array}{l}{\text{e}}^{-\frac{{\Vert x_i-x_j\Vert }^2}{\sigma }}{x}_j\in N_p\left(x_i\right)\\ 0其他\end{array}$ (11)

$L=D-W$， $D_{ii}={\displaystyle {\sum }_j{W}_{ij}}$，N_p(x_i)为数据点x_i的P个临近点， $\sigma$ 用来控制两个像元之间的相似程度，

一般取 $\sigma =1$。

根据梯度下降法，得到端元矩阵A和丰度矩阵S的更新法则。

$A\leftarrow A\frac{X{S}^{\text{T}}}{AS{S}^{\text{T}}+\alpha AQ}$ (12)

$S\leftarrow S\frac{A^{\text{T}}X+\mu SW}{A^{\text{T}}AS+\frac{1}{2}\lambda S^{-\frac{1}{2}}+\mu SD}$ (13)

3.2. 目标函数求解

目标函数初始化采用顶点成分分析法和全约束最小二乘算法相结合的方法(VCA-FCLS)。

其次，在迭代过程中，将观测矩阵X与端元矩阵A分别用Xt与At来替代，来保证丰度矩阵和为1约束。其定义如下：

$Xt=\left[\begin{array}{c}X\\ \delta 1_N^{\text{T}}\end{array}\right]$ (14)

$At=\left[\begin{array}{c}A\\ \delta 1_P^{\text{T}}\end{array}\right]$ (15)

其中 $\delta$ 为控制和为1约束效果的参数，一般将 $\delta$ 设置为25。

本次实验达到如下两个条件之一时，迭代停止：1) 当迭代次数到达最大；2) 当两次迭代结果之间的重构误差小于设定的阈值。如式(16)所示：

$\frac{\left|F_{K+1}-F_K\right|}{F_{K+1}}\le \tau$ (16)

$F_k$、 $F_{k+1}$ 为连续两次迭代的图像重构误差， $\tau$ 为设定的阈值。

3.3. 算法流程

DLGNMF算法具体步骤如表1所示。

Table 1. Detailed steps of DLGNMF algorithm unmixing

表1. DLGNMF算法具体步骤

4. 实验分析

本文采用仿真实验数据与真实高光谱数据来评价DLGNMF算法的解混性能，并与VCA-FCLS [9] 、MVC-NMF [10] 和GLNMF [10] 三种算法进行对比分析。

4.1. 评价指标

采用三个定量评价指标对解混性能行分析与评价 [9] [10]。

1) 光谱角距离(SAD)对提取的端元矩阵进行评价，其公式如(17)所示：

${\text{SA}}_i\text{D}={\cos }^{-1}\left(\frac{A_i^{\text{T}}{A}_i}{\Vert A_i\Vert \Vert {\hat{A}}_i\Vert }\right)$ (17)

其中 $A_i$， ${\hat{A}}_i$ 分别表示解混得到的端元光谱和真实端元光谱。

2) 均方根误差RMSE对仿真实验中反演的丰度矩阵进行评价，其公式如(18)所示：

$\text{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\Vert S_j-{\hat{S}}_j\Vert }^2}}$ (18)

其中 $S_j$， ${\hat{S}}_j$ 分别表示解混得到像元的丰度和真实图像中像元的丰度。

3) 由于本文的真实数据中没有丰度的先验知识，因此采用光谱重构误差(SRE)对真实高光谱图像的解混效果进行整体评价，其公式如(19)所示：

$\text{SRE}=\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\sqrt{\frac{1}{L}{\Vert X_j-{\hat{X}}_j\Vert }_2^2}}$ (19)

其中 $X_j$， ${\hat{X}}_j$ 分别表示重构观测矩阵的像元光谱和真实观测矩阵的像元光谱。

4.2. 仿真实验数据

从USGS光谱库选取5种纯净矿物作为端元成分，分别为Kaolinite#1 (高岭石#1)、Kaolinite#1 (高岭石#2)、Actinolite (阳起石)、Axinite (斧石)、Biotite (黑云母)，并将其记为P1、P2、P3、P4、P5。其光谱范围0.4~2.5 μm，波段数为224。其中高岭石#1，高岭石#2两种端元的光谱角距离SAD值为0.1299，说明两种端元的光谱曲线特征较为相似，以验证本文算法是否可以在“异物同谱”端元的情况下更好的提高解混精度。仿真图像大小为100 × 100，端元在仿真图像中的分布服从Matern高斯分布 [11] [12] [13] [14] [15]。本次仿真数据不包含纯像元，因此将每个端元对应的丰度值随机分布在0.05~0.8之间，最终得到无噪声的仿真高光谱图像。5种矿物光谱曲线和端元丰度分布如图2和图3所示。

Figure 2. Endmember spectral curves of five minerals

图2. 五种矿物端元光谱曲线

Figure 3. Abundance of five mineral endmember

图3. 五种矿物端元对应丰度

4.3. 参数的选择

采用控制单一变量的方法来确定最佳参数 $\alpha$、 $\lambda$、 $\mu$，即每次固定两个参数变换一个参数，当确定了最佳参数后，后面实验均采用本次最优参数值。

本次实验数据添加了SNR为30 dB的高斯白噪声。 $\lambda$ 是由图像自身的稀疏性决定 [9]，计算如式(20)所示。由仿真数据可得 $\lambda =0.37$。

$\lambda =\frac{1}{\sqrt{L}}{\displaystyle {\sum }_l\frac{\sqrt{N}-{\Vert x_l\Vert }_1-{\Vert x_l\Vert }_2}{\sqrt{N-1}}}$ (20)

首先设定 $\delta =25$、 $\sigma =1$、 $\lambda =0.37$、 $\mu =0.1$， $\alpha$ 从0.1变换到1，步长为0.1。其次固定 $\delta$、 $\sigma$、 $\lambda$ 不变， $\alpha =0.1$， $\mu$ 从0.1变换到1，步长为0.1。

在不同 $\alpha$ 和不同 $\mu$ 的情况下，DLGNMF算法所求得的SAD值和RMSE值如图4、图5所示。

Figure 4. Variation trend of SAD and RMSE with α

图4. SAD和RMSE随α的变换情况图

Figure 5. Variation trend of SAD and RMSE with μ

图5. SAD和RMSE随μ的变换情况图

通过图4，图5可知，本文算法在 $\alpha =0.1$、 $\mu =0.5$ 时，解混的效果最好，因此在后续的实验中，最优参数采用 $\delta =25$、 $\sigma =1$、 $\alpha =0.1$、 $\mu =0.5$， $\lambda$ 可根据公式(20)计算获得。其他对比算法的最优参数采用其论文中的最优参数设置。

此外，在流行正则约束中，构造像元的最近邻图时，需设定邻域窗口的大小。根据Cai等人的研究 [8]，算法在邻域窗口为5 × 5时效果最好。

4.4. 算法抗噪性检验

对模拟图像中添加信噪比为10 dB、20 dB、30 dB、40 dB、50 dB的高斯白噪声 [16]，以检验DLGNMF算法的抗噪性。对不同噪声下的图像进行解混，并与VCA-FCLS、MVC-NMF，GLNMF算法进行对比，以检验DLGNMF算法解混性能的优越性。

由于初始化的方法具有一定的随机性，取30次实验结果的平均值为最终解混结果。实验结果如表2、图6所示。

Table 2. Results of kaolinite#1 and kaolinite#2 under different algorithms

表2. 高岭石#1和高岭石#2在不同算法下的解混结果

Figure 6. Comparison of spectral curves of kaolinite#1 and kaolinite#2 obtained with real endmembers

图6. 高岭石#1、高岭石#2光谱曲线与真实端元的对比

从定量角度分析，表2中的数据可知，DLGNMF算法得到的高岭石#1、高岭石#2与真实端元之间的光谱角距离SAD值都是最小的。

从定性角度分析，图6中可知，DLGNMF得到的高岭石#1、高岭石#2端元光谱与真实端元的光谱曲线基本一致。说明DLGNNF算法可以更好的提高“异物同谱”端元存在情况下的解混精度，较为精确的提取出具有相似光谱曲线的端元。

Figure 7. SAD curves and RMSE curves of different algorithms under different SNR conditions

图7. 不同算法在不同信噪比条件下SAD曲线图与RMSE曲线图

从图7中可知，相比于其他算法，在不同信噪比的条件下，DLGNMF算法的解混精度最好，性能始终最优，说明该算法具有良好的抗噪性。

4.5. 图像大小对解混效果的影响

为测试不同图像大小对解混精度的影响，对图像添加30 dB的高斯白噪声，并将图像大小分别设置为60 × 60、80 × 80、100 × 100、120 × 120、140 × 140，结果如图8所示 [17]。

Figure 8. SAD curves and RMSE curves of different algorithms under different image sizes

图8. 不同算法在不同图像大小条件下SAD曲线图与RMSE曲线图

从图8中可以看出，相比于其他算法，DLGNMF算法在面对不同图像大小情况下，都展现出最优的解混精度。同时，随着图像的增大，SAD曲线与RMSE曲线趋势均下降，说明像元数目的增加对端元矩阵的提取与丰度矩阵的反演均有改善，可以提高解混的精度。

5. 真实高光谱数据实验验证

本次真实数据采用AVIRIS成像光谱仪拍摄的美国内达华州Cuprite矿区数据，是混合像元分解中运用较多的真实髙光谱数据 [18] [19] [20] [21]，如图9所示。该图像包含了224个波段，图像大小为250 × 190，光谱分辨率为10 nm，波长为0.38~2.5 μm。根据先验知识可知，该高光谱图像中主要包含9种端元成分，即Kaolinite1 (高岭石#1)、Kaolinite2 (高岭石#2)、Alunite (明矾石)、Dumortierite (蓝线石)、Nontronite (绿脱石)、Pyrope (镁铝榴石)、Montmorillonite (蒙脱石)，Chalcedony (玉髓)、Muscovite (白云母) [14]。考虑到有些波段的信噪比很低以及水蒸气吸收的影响(如波段1~3，105~115，150~170) [15]。在解混之前将这些波段剔除，留下189个波段进行数据的处理，参数设置与仿真实验相同。

Figure 9. AVIRIS hyperspectral data

图9. AVIRIS高光谱数据

图10表示DLGNMF算法所得到的部分端元的光谱。图11表示DLGNMF算法所得到9种不同端元的丰度图。做进一步的定量分析，表3所展示的是不同算法提取出的端元与真实端元光谱之间的比较，以及重构高光谱图像与真实高光谱图像之间的比较，分别用SAD与SRE表示。

Figure 10. Spectral curves corresponding to some endmembers. (a) Kaolinite#1; (b) Kaolinite#2; (c) Nontronite

图10. 部分端元对应的光谱曲线。(a) 高岭石#1；(b) 高岭石#2；(c) 绿脱石

Figure 11. Abundance corresponding to different endmembers. (a) Alunite; (b) Andradite; (c) Buddingtonite; (d) Dumortierite; (e) Kaolinite#1; (f) Kaolinite#2; (g) Montmorillonite; (h) Nontronite; (i) Chalcedony

图11. 不同端元所对应的丰度。(a) 明矾石；(b) 钙铁榴石；(c) 水铵长石；(d) 高岭石#1；(e) 高岭石#2；(f) 蓝线石；(g) 绿脱石；(h) 蒙脱石；(i) 玉髓

从图10可知，DLGNMF算法得到的高岭石#1、高岭石#2光谱曲线与真实地物的光谱曲线基本一致，较好的将相似端元光谱解混出来，DLGNNF算法可以更好的提高“异物同谱”端元存在情况下的解混精度，较为精确的提取出具有相似光谱曲线的端元。

从表3中可以看出，相比于其他算法，DLGNMF算法提取的端元光谱曲线，其SAD值大部分都是最小的，SRE值也是最小。说明在相同条件下，该算法的解混精度较高，性能较好，优于其他同类算法。

Table 3. SAD and values obtained by different algorithms for this region

表3. 不同算法对该地区解混得到的SAD值与SRE值

6. 结论

针对非负矩阵解混模型容易陷入局部极小值、受噪声影响较大以及混合像元中存在“异物同谱”端元的情况下解混精度较低的问题。本文将空间信息与光谱信息相结合，在非负矩阵解混模型的基础上添加单形体最小距离的端元约束与稀疏图正则的丰度约束，提出了基于最小距离和稀疏图正则约束的非负矩阵解混算法(DLGNMF)。通过仿真数据与真实数据进行算法的性能测试，并与VCA-FCLS、GLNMF、MVC-MNF三种算法进行对比。实验结果表明：

相比于其他同类算法，DLGNMF算法弥补了NMF算法解混过程中高光谱数据空间几何特性的缺失，能够较为精确的提取出具有“异物同谱”的端元成分。在SAD、RMSE、SRE方面均获得最小值，表明DLGNMF算法解混精度较高，具有更好的解混性能。

DLGNMF算法同时兼顾了图像的全局稀疏性与局部相似性，使得到的丰度矩阵更加符合真实端元分布。即使在处理噪声较大的高光谱数据时，解混精度更高，抗噪性能更好，结果更加稳定。

需要指出，添加流行正则约束时，矩阵W、D、L，都是N × N大小的矩阵(N为像元个数)，使得迭代过程中矩阵占据过大的电脑内存，严重影响算法的运算效率。因此如何在提高解混精度的同时提升运算效率，是后续重点的研究内容。

NOTES

^*通讯作者。