Cauchy Nonnegative Matrix Factorization for Hyperspectral Unmixing Based on Graph Laplacian Regularization

Shanxue Chen; Shaohua Xu

doi:10.3788/LOP232437

Key worlds hyperspectral unmixing; nonnegative matrix factorization; Cauchy loss function; graph Laplacian; reweight sparse

1　引言

由成像光谱仪^［1］采集到的具有数百或数千波段的高光谱图像（HSI）已经成为多个领域的研究热点。凭借丰富的光谱信息，高光谱图像已经广泛应用于农业、食品质量控制、矿物勘探、生物医学等诸多热门领域^［2-4］。然而，受仪器性能的限制、时空复杂度的影响，较低的分辨率会使单个像元可能含有多种成分，导致混合像元的产生。为了解决这个问题，通常使用高光谱解混来分解混合像元，将高光谱图像的每个像元分解为一组端成员及其相应的丰度组合。

近年来，基于线性混合模型（LMM）^［5］和欧氏距离的标准非负矩阵分解（NMF）^［6］模型备受广大研究者的关注，它是一种盲源分解技术。由于目标函数的非凸性，大量的局部最小值导致目标式陷入局部最优解。因此，常常在NMF模型中通过端元和丰度的先验信息添加一些合适的约束以缩减解空间。

关于端元的约束：从几何角度出发，Miao等^［7］提出最小体积约束的NMF（MVCNMF），它将最小体积约束纳入分解过程中，从而减小估计数据与真实端元之间的散度。Pompili等^［8-9］将正交约束加入NMF中，在保证端元间相互独立的同时也减小了噪声对解混结果的影响。Jia等^［10］将分段光滑约束加入NMF中以防止光谱数据出现不连续。另一方面，丰度具有更加丰富的先验信息，对其添加约束会起到事半功倍的效果。Qian等^［11］、Chen等^［12］分别利用L_1/2正则项、双空间加权稀疏约束项来促进丰度的稀疏并取得了较好的效果。Ince^［13］在保证丰度稀疏的基础上充分挖掘HSI数据中潜在的流形结构，提出了双空间的基于图拉普拉斯正则化的稀疏解混算法并取得一定的效果。然而，传统的NMF对噪声和异常像元非常敏感。由于高斯白噪声和椒盐噪声的存在，基于欧氏距离的传统NMF很容易丢失精度，造成解混的失败。针对以上问题，陈善学等^［14］将L₂₁范数引入标准的NMF模型中并且使用双重加权稀疏约束和子空间结构正则化对丰度进行约束，具有较好的鲁棒性。吴新浪等^［15］采用Cauchy非负张量分解代替传统NMF方法来降低噪声影响，同时又保留了HSI的空间结构信息，取得较好的解混效果。

为了更加有效地利用空间信息，提高鲁棒性和解混性能，本文提出一种新的基于图拉普拉斯正则化的柯西非负矩阵分解（CNMF-GLR）算法。首先，采用基于柯西损失函数模型以获得良好的鲁棒性。其次，为提高解混性能防止丰度的内在结构遭到破坏，将重加权稀疏约束与图拉普拉斯约束加入模型中。由于图拉普拉斯算法在HSI上对于给定的像元只是将欧氏空间上相邻的像元相连^［16］，因此相邻点的选择便成了关键。于是，使用一种局部邻域加权方法，通过矩形窗确定目标像元的邻域，构建图拉普拉斯正则（GLR）对局部邻域的流形结构进行学习。最后通过实验验证所提算法的有效性。

2　基本原理

2.1　线性混合模型与非负矩阵分解

通常情况下，观测到的HSI通常是一个矩阵 $Y \in R^{L \times N}$ ，在LMM下第k个像元 $y_{k}$ 可以表示为端元光谱和其对应比例相乘的形式，即

y_{k} = \sum_{i = 1} e_{i} a_{i} + n,

（1）

式中： $e_{i}$ 为第i个端元光谱； $a_{i}$ 为对应的丰度； $n$ 为噪声。忽略噪声的影响，将式（1）写成矩阵的形式，则有

Y \approx E A

，（2）

式中： $E \in R^{L \times K}$ ，是端元矩阵； $A \in R^{K \times N}$ ，为丰度矩阵；N、L、K分别表示为HSI像元的总数目、波段数量和端元个数。为了衡量近似差异，通过欧氏距离来计算Y与EA的差距，考虑非负性和丰度和唯一性可建立NMF模型：

\underset{E, A}{a r g m i n} {‖Y - E A‖}_{F}^{2} s . t . E \geq 0, A \geq 0, 1_{k}^{T} A = 1_{N}^{T} 。

（3）

Lee等^［17］通过梯度下降法给出了式（3）的乘法更新法则：

E \leftarrow E . * (Y A^{T}) . / (E A A^{T})

，（4）

A \leftarrow A . * (E^{T} Y) . / (E^{T} E A)

，（5）

式中：.*表示矩阵对应元素相乘；./表示对应元素相除。

2.2　Cauchy损失函数

柯西损失函数（CLF）^［18］的表达式为

ϕ (x) = l n [1 + {(\frac{x}{r})}^{2}]

，（6）

式中：r为可调常数。基于式（6）的解混模型虽然对鲁棒性有一定的提升，但是对离群值并没有那么敏感。极端情况下，无限大的噪声会产生难以估计的能量。针对这一问题，文献［19］对该模型进行改进并提出基于截断柯西损失函数的模型，改进后的柯西损失函数为

f (x) = \{\begin{matrix} l n [1 + {(x / r)}^{2}], |x| \leq c \\ l n [1 + {(c / r)}^{2}], |x| > c \end{matrix}

。（7）

从式（7）不难看出，当 $x$ 大于阈值c时会用c来代替损失函数的值从而达到抑制和平滑异常值的目的，提高了NMF的稳定性。图1给出了传统损失函数与截断式损失函数的比较。

Figure 1.Comparison of different loss function models

Download full size

View all figures

图1中横轴 $x$ 为损失函数的自变量，纵轴为损失值^［19］。从图中可以较直观观察到，截断式柯西损失对异常值具有较好的鲁棒作用。基于截断式柯西损失函数模型可以描述为

\underset{E, A}{a r g m i n} \frac{1}{2} \sum_{i = 1}^{L} \sum_{j = 1}^{N} f [{(Y - E A)}_{i j}] 。

（8）

3　基于图拉普拉斯正则化的柯西非负矩阵分解模型

3.1　CNMF-GLR模型

为提高NMF的鲁棒性，抑制异常像元对分解的影响，采用柯西损失函数模型来代替传统NMF模型，有

J (E, A) = \frac{1}{2} \sum_{i = 1}^{L} \sum_{j = 1}^{N} f [{(Y - E A)}_{i j}] 。

（9）

大多数情况下，HSI中的像元不会含有所有端元成分。因此，丰度矩阵是稀疏的。为缩减NMF的解空间，提高解混精度，引入重加权稀疏约束^［10］：

\underset{}{m i n} {‖Q \otimes A‖}_{_{1}} s . t . Y = E A,

（10）

式中： $\otimes$ 是哈达玛积，表示两个矩阵对应元素相乘。自适应权重矩阵 $Q$ 定义为

Q^{t + 1} = 1 . / (A^{t} + C_{e p s}),

（11）

式中， $C_{e p s}$ 为极小正常数矩阵，用来防止分母为0。

将重加权稀疏约束项加入CNMF模型中，则有

J (E, A) = \frac{1}{2} \sum_{i = 1}^{L} \sum_{j = 1}^{N} f [{(Y - E A)}_{i j}] + λ_{1} {‖Q \otimes A‖}_{1} 。

（12）

为了挖掘HSI内在流形结构，可以将HSI数据看成是位于一个流形空间中。将HSI中的每一个像元视为一个顶点，N个像元可以构成一个N×N的邻接矩阵S。注意到，一个完整图的边数与节点数量呈二次关系^［16］，因此HSI的NMF可能会变得非常复杂。为了解决这个问题，考虑HSI邻近的数据点在新的数据空间中也是相近的，假设HSI数据位于高维流形空间中的低维曲面上，参考文献［20］可以使用图拉普拉斯映射的方法设置权重。对于给定的N个数据点的数据集 $Y = {y_{1}, y_{2}, \dots, y_{N}}$ 构建邻接矩阵 $S \in R^{N \times N}$ ，有

S_{i j} = \{\begin{array}{l} e x p (\frac{- {‖y_{i} - y_{j}‖}^{2}}{2 θ^{2}}), & i i s a d j a c e n t t o j \\ 0, & o t h e r w i s e \end{array}

，（13）

式中： $θ$ 为权衡参数，用于控制核值的大小。将HSI矩阵Y映射为丰度矩阵A， $y_{i}$ 与 $y_{j}$ 越相似， $S_{i j}$ 越接近1， $a_{i}$ 与 $a_{j}$ 也越相近。因此图拉普拉斯正则化可以表述为

m i n t r (A L A^{T}) = \frac{1}{2} \sum_{i, j} S_{i j} {‖A_{i} - A_{j}‖}_{2}^{2}

，（14）

式中：tr（·）表示矩阵的迹； $L = D - S$ ， $D$ 为对角矩阵， $D_{i i} = \sum_{j = 1}^{N} S_{i j}$ 。

由于使用Cauchy损失模型可能会破坏像元之间的相关性，而像元之间的相关性在流形空间中可以映射为丰度之间的联系。因此，为了保持丰度矩阵内在相关性不被破坏，提高解混性能，将图拉普拉斯正则化加入模型中，提出一种新的CNMF-GLR模型：

\begin{array}{l} m i n J (E, A) = \frac{1}{2} \sum_{i = 1}^{L} \sum_{j = 1}^{N} f [{(Y - E A)}_{i j}] + λ_{1} {‖Q \otimes A‖}_{1} + \\ λ_{2} t r (A L A^{T}) s . t . E \geq 0, A \geq 0, 1_{k}^{T} A = 1_{N}^{T} \end{array}

。（15）

式中：正则化参数 $λ_{1} > 0$ 、 $λ_{2} > 0$ 。

3.2　模型优化与迭代公式推导

由于式（15）是非线性、非凸的，很难直接求解。可以使用半二次技术^［21］来解决优化问题（15）。假设 $g (x) = - f (x)$ ，引入 $φ (X)$ 作为g（x）的共轭函数，其中，X是相应的共轭变量，采用坐标下降法^［22］将式（15）优化：

\begin{array}{l} m i n J (E, A) = \sum_{i j} X_{i j} \frac{{(Y - E A)}_{i j}^{2}}{2 r^{2}} + λ_{1} {‖Q \otimes A‖}_{1} + \\ λ_{2} t r (A L A^{T}) s . t . E \geq 0, A \geq 0, 1_{k}^{T} A = 1_{N}^{T} \end{array}

。（16）

由文献［19］可知，当E、A固定时，式（16）中X的最优解为

X_{i j} = \{\begin{array}{l} \frac{1}{1 + \frac{{(Y - E A)}_{i j}^{2}}{r^{2}}}, {|Y - E A|}_{i j} \leq r \sqrt[]{c} \\ 0, {|Y - E A|}_{i j} > r \sqrt[]{c} \end{array},

（17）

可以使用乘性迭代规则来解决优化式（16）。

对于目标函数（16）引入拉格朗日乘子 $Γ$ ：

\begin{array}{l} J (A) = \sum_{i j} {X_{f}}_{i j} \frac{{(Y_{f} - E_{f} A)}_{i j}^{2}}{2 r^{2}} + λ_{1} {‖Q \otimes A‖}_{1} + \\ λ_{2} t r (A L A^{T}) + T r (Γ A^{T}) \end{array}

（18）

式中： $Y_{f}$ 、 $A_{f}$ 和 $X_{f}$ 是增广矩阵：

\{\begin{matrix} \begin{matrix} Y_{f} = [\begin{matrix} Y \\ δ 1_{N}^{T} \end{matrix}] \\ A_{f} = [\begin{matrix} A \\ δ 1_{N}^{T} \end{matrix}] \end{matrix} \\ X_{f} = [\begin{matrix} X \\ δ 1_{N}^{T} \end{matrix}] \end{matrix},

（19）

式中： $δ$ 为正常数，是用于调节丰度和为1的参数^［23］。对矩阵 $E$ 求偏导，在库恩-塔尔条件（KKT）下，有

\begin{array}{l} \nabla_{A} J (A) = - E^{T} (X_{f} \otimes Y_{f}) + E_{f}^{T} (X_{f} \otimes E_{f} A) + λ_{1} Q - \\ λ_{2} A S + λ_{2} A D + Γ = 0 \end{array}

，（20）

A \otimes Γ = 0,

（21）

对式（20）两侧同时 $\otimes A$ 并将式（21）代入可得 $A$ 的更新规则：

A \leftarrow A \otimes \{[E_{f}^{T} (X_{f} \otimes Y_{f}) + λ_{2} A S] . / [E_{f}^{T} (X_{f} \otimes E_{f} A) + λ_{1} Q + λ_{2} A D]\},

（22）

对于式（20），

- λ_{2} A S + λ_{2} A D,

（23）

它是由式（14）对A求偏导而来的。通过式（13），可以观察到矩阵 $S, D \in R^{N \times N}$ 的求解复杂度较高且与像元邻域相关，在求解式（23）时由于 $A \in R^{K \times N}$ 使得时间复杂度较高且属于高数量级 $O (K \times N^{2}) \approx O (N^{2})$ 。可以通过局部邻域加权^［24］的方法确定局部邻域以减少时间复杂度。具体做法如下：

对于给定的像元 $y_{i}$ ，图拉普拉斯映射只考虑 $y_{i}$ 邻近像元对其的贡献。首先，将HSI数据矩阵 $Y \in R^{L \times N}$ 的N个像元排列在一个 $m \times n$ 的网格中。像元 $y_{i}$ 与它网格的位置 $(c, l)$ 关系可以表示为

i = (c - 1) m + l,

（24）

i = (l - 1) n + c 。

（25）

其次，将丰度矩阵 $A \in R^{K \times N}$ 以同样的方式排列在网格中。 $y_{i}$ 局部邻域内的全部像元对它的贡献 $w_{i}$ 可以表示为

w_{i} = \sum_{j \in N (j)} S_{i j}, i = 1,2, 3, \dots, N 。

（26）

因此，矩阵 $W \in R^{K \times N}$ 可以构建为

W = {[\begin{matrix} w_{1} & w_{2} & \dots & w_{N} \\ w_{1} & w_{2} & \dots & w_{N} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ w_{1} & w_{2} & \dots & w_{N} \end{matrix}]}_{K \times N} 。

（27）

对于邻域的选取，参考文献［13］、［22］和［24］，可以选择矩形结构作为局部邻域。每一个中心像元的丰度 $A_{i}$ 如图2所示，可以通过邻域像元的丰度对 $A_{i}$ 的贡献去评估 $A_{i}$ 与周围像元丰度的关系。不难得到：

Figure 2.Distribution of pixel and abundance in local neighborhood

Download full size

View all figures

A S = A_{W}

，（28）

A D = W \otimes A

，（29）

式中： $A_{W} \in R^{K \times N}$ ， $A_{W}$ 的第i列可以表示为 ${\tilde{a}}_{i} = \sum_{j \in N (i)} S_{i j} {\tilde{a}}_{j}$ 。矩阵S与W可在算法初始化时求解并作为结果常量直接参与迭代过程。所以迭代过程中求解 $A_{W}$ 矩阵只与邻域中像元数量有关。由于邻域像元数量为常数级，因此可以得到式（28）的时间复杂度为

\begin{array}{l} O (G \times m \times n) = O (G \times \sqrt[]{N} \times \sqrt[]{N}) = \\ O (G \times N) \approx O (N) \end{array}

，（30）

式中：G为一常数，表示邻域像元的数目。由于式（29）的时间复杂度为 $O (K \times N) \approx O (K)$

所以，式（28）和式（29）的时间复杂度为 $O (N)$ ，优于式（23）。因此，将式（28），（29）代入式（22）中可得A的更优更新规则：

A \leftarrow A \otimes \{[E_{f}^{T} (X_{f} \otimes Y_{f}) + λ_{2} A_{w}] . / [E_{f}^{T} (X_{f} \otimes E_{f} A) + λ_{1} Q + λ_{2} W \otimes A]\}

。（31）

对端元矩阵E的估计，可以对目标函数（16）引入拉格朗日乘子 $Ψ$ ，有

m i n J (A) = \sum_{i j} X_{i j} \frac{{(Y - E A)}_{i j}^{2}}{2 r^{2}} + T r (Ψ A^{T})

，（32）

在KKT条件下，有

\nabla_{E} J (E) = - (X \otimes Y) A^{T} + (X \otimes E A) A^{T} + Ψ = 0

，（33）

E \otimes Ψ = 0

，（34）

对式（33）两侧同时 $\otimes E$ 可获得端元 $E$ 更新公式：

E \leftarrow E \otimes [(X \otimes Y) A^{T}] . / [(X \otimes E A) A^{T}]

。（35）

3.3　算法实现

在优化问题中，初始化方式的选择非常关键，合适的初始化方法可以有效缩减解空间，提升收敛速率。对于HSI解混问题通常有两种初始化方法：随机初始化和顶点成分分析-完全约束最小二乘（VCA-FCLS）初始化^［24］。使用VCA初始化端元矩阵E，FCLS初始化丰度矩阵A可以获得更高的精度和解混效率，所有实验均采用VCA-FCLS初始化的方式。

为获得最优的局部流形学习效果，通过实验选择合适大小的矩形窗结构作为图拉普拉斯的局部邻域进行学习。

对于式（19），参考文献［25］将 $δ$ 设置为18。CNMF-GLR的迭代终止条件有两个：1）最大迭代次数；2）两次迭代重构图像之间的欧氏距离小于某一阈值。

d = {‖W^{t + 1} H^{t + 1} - W^{t} H^{t}‖}_{F}^{2}

。（36）

在所有实验中，最大迭代次数设置为3000，阈值设置为 $10^{- 4}$ 且连续10次小于该阈值则认为模型收敛，停止更新。具体流程如图3所示。

Figure 3.Flow of CNMF-GLR algorithm

Download full size

View all figures

4　实验结果与分析

为评估CNMF-GLR算法的有效性，在模拟数据集和真实数据集上进行实验，并与MVCNMF^［7］、CauchyNMF^［18］、L₂₁NMF-SSR^［14］、DSGLSU^［13］和SSCNMF^［25］算法进行对比。由于初始值对实验结果影响较大，采用同一初始化方法在不同实验下也会产生差异较大的初始化结果，难以有效评估。针对此问题，所有对比实验均采用相同的初始化矩阵E、A，通过不同模型的实验结果来评估算法的性能。通常情况下，使用光谱角距离（DD_SAD）来计算估计端元 $\hat{E}$ 与真实端元 $E$ 之间的相似程度，采用均方根误差（E_RMSE）来计算重构图像 $\hat{Y}$ 与真实图像 $Y$ 的相似程度。

{D_{S A D}}_{(e, \hat{e})} = a r c c o s (\frac{e^{T} \hat{e}}{{‖e‖}_{2} {‖\hat{e}‖}_{2}})

，（37）

{E_{R M S E}}_{(Y, \hat{Y})} = \frac{1}{N} \sum_{i = 1}^{N} \sqrt[]{\frac{1}{L} {‖y_{i} - {\hat{y}}_{i}‖}_{2}^{2}}

，（38）

式中： $e 、 \hat{e}$ 为端元向量； $y 、 \hat{y}$ 为像元向量。不难看出，D_SAD和E_RMSE的值越小说明解混精度更高，以此来评估算法性能。

所有实验均在相同环境下进行：使用的电脑运行内存为8 GB，处理器为i5-8500，主频率为3 GHz，使用MATLAB R2016a仿真平台，在Windows 10环境下运行。Jasper Ridge是一种广泛应用于HSI解混实验的数据集。原始数据集中记录了380~2500 nm中的224个通道，含有512×614个像元。由于该HSI数据过于庞大，难以处理，因此，只考虑原始数据集中100×100的子图像。为了避免大气和水蒸气影响，去掉了1~3、108~112、154~166和220~224波段，保留198个波段。数据集中同时也包含该HSI中4种端元光谱，它们分别为植被、水体、土壤和道路。相似地，Urban数据集也是一种广泛使用的HSI数据集，经过预处理，移除了1~4、76、87、10~111、136~153和198~210波段，保留162个波段，使用307×307个像元的子图像。该数据集中主要包含4种地物：沥青、草地、树木和屋顶。图4展示了Jasper Ridge数据集中4种地物的光谱图。其中，横坐标表示的是光谱通道，纵坐标表示的是归一化后的光谱强度。

Figure 4.Spectra of 4 ground objects in the Jasper Ridge dataset

Download full size

View all figures

4.1　模拟数据实验

使用MATLAB高光谱合成工具箱^［14］生成与真实场景图像相似的模拟图像。使用Jasper Ridge数据集中真实端元光谱作为实验光谱，根据球面高斯场生成相应的丰度，按照每个像元中各个端元的丰度将端元光谱混合，即可得到无噪声的模拟图像。由于真实HSI图像中普遍存在高斯白噪声和一些异常点，因此在模拟图像中添加了不同等级的高斯白噪声和椒盐噪声以验证所提算法的有效性。

4.1.1　参数设置

图5为所提算法在不同大小矩形窗下的表现。由于邻域像元数量为5×5与7×7矩形窗表现相近，考虑到计算复杂度的问题，后续实验邻域均选择5×5的矩形窗。同时，在所提模型中，有两个关键参数 $λ_{1}$ 和 $λ_{2}$ ， $λ_{1}$ 用来控制重加权稀疏约束项， $λ_{2}$ 控制拉普拉斯正则项。本小节主要研究这两个参数在不同取值下对算法性能的影响。为了更加接近真实数据在无噪声的模拟图像中加入信噪比（SNR）^［14］为25的高斯白噪声，加入2%的异常像元（等级为0.02的椒盐噪声）。 $λ_{1}$ 、 $λ_{2}$ 分别以｛0.0005，0.001，0.005，0.001，0.002，0.003｝为间隔，考虑它们的所有排列组合，每种组合进行10次实验取均值。

Figure 5.Comparison of average D_SAD and E_RMSE in different neighborhoods

Download full size

View all figures

图6展示了不同 $λ_{1}$ 和 $λ_{2}$ 在模拟数据集上的表现，从图6可以看出，D_SAD和E_RMSE的变化趋势具有一致性。当 $λ_{1}$ 为0.01， $λ_{2}$ 为0.2时，D_SAD和E_RMSE可以得到相对较优的结果。因此，所有实验 $λ_{1}$ 均设置为0.01， $λ_{2}$ 设置为0.2。

Figure 6.Average D_SAD and E_RMSE of CNMF-GLR under different $λ_{1}$ and $λ_{2}$ values

Download full size

View all figures

4.1.2　性能分析

为探究所提算法的鲁棒性和有效性，实验分为两组，其中，一组是在无噪声的模拟数据集中分别加入信噪比等级为10 dB、15 dB、20 dB、25 dB的高斯白噪声。表1和表2展示了不同SNR等级下各算法的D_SAD和E_RMSE的表现。从表1和表2中的数据可以看出，在不同等级高斯白噪声的影响下，CNMF-GLR算法相比于其他算法拥有相对较低的D_SAD值和E_RMSE值，说明所提算法具有不错的抗干扰能力。与无约束的CauchyNMF模型和使用稀疏约束的SSCNMF相比，利用流形空间结构的CNMF-GLR具有更优的表现，证明了图拉普拉斯正则化保持丰度流形结构的有效性。然而，随着SNR逐渐增大，无约束的CauchyNMF在低噪声的环境下精度提升反而不如经典算法MVCNMF，这说明Cauchy模型虽然具有较好的鲁棒性，但也可能破坏HSI数据的内在结构，导致解混性能降低。而所提算法使用图拉普拉斯正则化维护了HSI数据内部结构，仍然具有较强的性能，说明所提算法也适用于低噪声的环境。

Table 1. Comparison of average D_SAD of various algorithms under the influence of Gaussian white noise

View table

View all Tables

Table 1. Comparison of average D_SAD of various algorithms under the influence of Gaussian white noise

SNR /dB	MVCNMF /arb.units	CauchyNM-F /arb.units	L21NMF-SSR /arb.units	DSGLSU /arb.units	SSCNMF /arb.units	CNMF-GLR /arb.units
10	0.35536	0.30956	0.22786	0.24786	0.22897	0.20864
15	0.22466	0.20687	0.17657	0.19756	0.18864	0.17056
20	0.15284	0.13786	0.09764	0.10258	0.10189	0.09540
25	0.07867	0.07246	0.04556	0.07045	0.05874	0.04626
30	0.01897	0.03765	0.00761	0.01087	0.00996	0.00789

Table 2. E_RMSE comparison of various algorithms under the influence of Gaussian white noise

View table

View all Tables

Table 2. E_RMSE comparison of various algorithms under the influence of Gaussian white noise

SNR /dB	MVCNMF /arb.units	CauchyNM-F /arb.units	L21NMF-SSR/arb.units	DSGLSU/arb.units	SSCNMF / arb.units	CNMF-GLR / arb.units
10	0.087496	0.071213	0.069864	0.082456	0.070896	0.066233
15	0.057897	0.061012	0.047892	0.048964	0.052168	0.047223
20	0.048887	0.053396	0.034687	0.041268	0.039864	0.032899
25	0.028794	0.036765	0.013786	0.020489	0.018964	0.014896
30	0.012560	0.011587	0.006987	0.009645	0.008650	0.007231

在实际的高光谱图像中，往往会附带一些异常点。通常这些异常点是由图像传输中尖锐和干扰形成的。因此，另一组实验在无噪声的模拟数据集中分别加入密度等级为0.02、0.04、0.06、0.08、0.1的椒盐噪声。分析表3和表4中的数据可知，MV-CNMF受椒盐噪声影响较大。而采用柯西损失函数的算法SSCNMF与CNMF-GLR在面对受椒盐噪声污染的HSI时，拥有较低的D_SAD与E_RMSE值，这说明柯西模型在此环境下具有较为明显的优势。而CNMF-GLR模型表现出最低的数值，说明所提模型具有优越的鲁棒性和更优的性能，随着椒盐噪声的等级逐渐增大，CNMF-GLR的优势愈加明显。

Table 3. Comparison of average D_SAD of various algorithms under the influence of salt and pepper noise

View table

View all Tables

Table 3. Comparison of average D_SAD of various algorithms under the influence of salt and pepper noise

Grade/ arb.units	MVCNMF/ arb.units	CauchyNM-F/ arb.units	L₂₁NMF-SSR/ arb.units	DSGLSU/ arb.units	SSCNMF/ arb.units	CNMF-GLR/ arb.units
0.02	0.15814	0.12948	0.09764	0.10246	0.08979	0.08245
0.04	0.17687	0.13687	0.10788	0.13765	0.09864	0.08978
0.06	0.18094	0.15255	0.11868	0.15867	0.11736	0.09877
0.08	0.19924	0.16786	0.13254	0.17486	0.12387	0.11187
0.1	0.23948	0.16356	0.14786	0.18765	0.12278	0.11089

Table 4. E_RMSE comparison of various algorithms under the influence of salt and pepper noise

View table

View all Tables

Table 4. E_RMSE comparison of various algorithms under the influence of salt and pepper noise

Grade /arb.units	MVCNMF /arb.units	CauchyNM-F /arb.units	L₂₁NMF-SSR /arb.units	DSGLSU /arb.units	SSCNMF / arb.units	CNMF-GLR / arb.units
0.02	0.092149	0.068497	0.056798	0.078465	0.061286	0.039765
0.04	0.146836	0.087658	0.062488	0.098654	0.062186	0.057865
0.06	0.220817	0.113244	0.128247	0.143284	0.087657	0.064692
0.08	0.301487	0.154658	0.133684	0.167687	0.092468	0.070258
0.1	0.387996	0.155765	0.138994	0.193245	0.102865	0.079864

4.2　真实数据实验

图7（a）和图7（b）分别是Jasper Ridge数据集、Urban数据集在波段［198，34，28］和波段［38，28，12］的RGB显示。

Figure 7.RGB display on the real data set

Download full size

View all figures

图8显示了CNMF-GLR在Jasper Ridge数据集上提取到的端元光谱与参考光谱的对比，不难看出，提取端元光谱与真实光谱相近，这说明所提算法提取端元光谱较为有效。表5展现了不同算法下的D_SAD数据与运行时间。从表中可以看出，虽然所提算法运行时间不占优势，但对4种地物端元得到的D_SAD数值最小，说明CNMF-GLR算法牺牲了部分时间提取的端元光谱更加接近真实端元。表6对不同算法下的E_RMSE值进行对比，CNMF-GLR拥有最低的E_RMSE值，说明所提算法重构的图像与真实图像更加接近。不同算法丰度图像如图9所示，可以观察出，CNMF-GLR算法提取到的丰度图像轮廓最为清晰，进一步说明了该算法的有效性。

Figure 8.Comparison between end element spectra extracted by CNMF-GLR and reference spectra

Download full size

View all figures

Table 5. Comparison of D_SAD values of different algorithms in the Jasper Ridge dataset

View table

View all Tables

Table 5. Comparison of D_SAD values of different algorithms in the Jasper Ridge dataset

End member	MVCNMF / arb.units	CauchyNM-F / arb.units	L21NMF-SSR / arb.units	DSGLSU / arb.units	SSCNMF / arb.units	CNMF-GLR/ arb.units
Runtime /min	8.76	10.24	12.64	13.25	14.37	13.75
water	0.13989	0.08765	0.06276	0.07265	0.07686	0.03213
vegetation	0.12265	0.12023	0.11082	0.12013	0.11265	0.06025
soil	0.17894	0.17214	0.15786	0.16896	0.15876	0.11611
road	0.16379	0.09065	0.06686	0.09685	0.09241	0.05312
mean	0.15131	0.11766	0.09957	0.11464	0.11017	0.06541

Table 6. Comparison of average E_RMSE values of different algorithms in the Jasper Ridge dataset

View table
View all Tables
Table 6. Comparison of average E_RMSE values of different algorithms in the Jasper Ridge dataset

MVCNMF /
arb.units
CauchyNM-F /arb.units
L₂₁NMF-SSR /
arb.units
DSGLSU /
arb.units
SSCNMF /
arb.units
CNMF-GLR /
arb.units
0.067262 0.047274 0.031131 0.042862 0.038652 0.019101

Figure 9.Abundance image comparison of different algorithms in the Jasper Ridge dataset

Download full size

View all figures

另一组实验是在Urban数据集上进行的。表7和表8展现了不同算法的D_SAD值和E_RMSE值在该数据集下的表现。在表7中，由于该数据集中的数据多于Jasper Ridge数据集，因此各算法运行时间明显变长。通过表7对比D_SAD数据可以看出，对于草地和树木，CNMF-GLR具有最低的数值，说明对这两种端元提取到的光谱与真实光谱更加接近，总体来看表现更优。图10可以更为直观地反映这一点。在表8中使用CNMF-GLR重构的HSI图像拥有最低的E_RMSE值，说明更加接近真实图像。从图11丰度图像上可以观察到，CNMF-GLR提取到的各地物轮廓较为清晰，说明了所提算法的有效性。

Table 7. Comparison of D_SAD values of different algorithms in Urban dataset

View table

View all Tables

Table 7. Comparison of D_SAD values of different algorithms in Urban dataset

End member	MVCNMF / arb.units	CauchyNM-F /arb.units	L₂₁NMF-SSR / arb.units	DSGLSU / arb.units	SSCNMF / arb.units	CNMF-GLR / arb.units
Runtime /min	12.15	15.76	16.88	16.74	18.24	17.68
asphalt	0.34587	0.30221	0.22047	0.24521	0.10275	0.11346
grass	0.32786	0.31865	0.25738	0.30125	0.05860	0.03061
tree	0.20166	0.21636	0.13104	0.20106	0.03998	0.03284
roof	0.22765	0.20168	0.14107	0.19201	0.14687	0.14558
mean	0.27576	0.25972	0.18749	0.23489	0.08705	0.08062

Table 8. Comparison of average E_RMSE values of different algorithms in Urban dataset

View table
View all Tables
Table 8. Comparison of average E_RMSE values of different algorithms in Urban dataset

MVCNMF /
arb.units
CauchyNM-F /arb.units
L₂₁NMF-SSR /
arb.units
DSGLSU /
arb.units
SSCNMF /
arb.units
CNMF-GLR /
arb.units
0.28365 0.26542 0.132851 0.18652 0.082145 0.064547

Figure 10.Comparison between end element spectra extracted by CNMF-GLR and reference spectra

Download full size

View all figures

Figure 11.Abundance image comparison of different algorithms in Urban dataset

Download full size

View all figures

5　结论

提出一种新的高光谱解混算法CNMF-GLR。实验结果表明，只采用柯西损失函数的解混模型，虽然会获得较好的鲁棒性，但是会破坏数据的内在结构，导致精度丢失，在低噪声的环境下反而表现不佳。通过施加图拉普拉斯约束与重加权稀疏约束，可以有效避免这个问题。对于图拉普拉斯流形正则项，使用一种新的基于局部邻域加权的方法，在减小计算复杂度的同时也具备更加明确的物理意义。在模拟数据集、Jasper Ridge和Urban真实数据集上的实验结果表明，CNMF-GLR算法用于高光谱解混时具有更优的性能。然而，所提算法依然有一些不足：1）CNMF-GLR仅利用了丰度的先验信息，未考虑对端元施加约束；2）控制约束的正则化参数 $λ_{1}$ 和 $λ_{2}$ 无法自适应选择，在不同场景下合适的正则化参数可能会带来更大的收益。

Category: Remote Sensing and Sensors

Received: Nov. 6, 2023

Accepted: Nov. 27, 2023

Published Online: Jul. 8, 2024

The Author Email: Shaohua Xu (535849135@qq.com)

DOI:10.3788/LOP232437

CSTR:32186.14.LOP232437

Table 1. Comparison of average DSAD of various algorithms under the influence of Gaussian white noise

Table 1. Comparison of average DSAD of various algorithms under the influence of Gaussian white noise

Table 2. ERMSE comparison of various algorithms under the influence of Gaussian white noise

Table 2. ERMSE comparison of various algorithms under the influence of Gaussian white noise

Table 3. Comparison of average DSAD of various algorithms under the influence of salt and pepper noise

Table 3. Comparison of average DSAD of various algorithms under the influence of salt and pepper noise

Table 4. ERMSE comparison of various algorithms under the influence of salt and pepper noise

Table 4. ERMSE comparison of various algorithms under the influence of salt and pepper noise

Table 5. Comparison of DSAD values of different algorithms in the Jasper Ridge dataset

Table 5. Comparison of DSAD values of different algorithms in the Jasper Ridge dataset

Table 6. Comparison of average ERMSE values of different algorithms in the Jasper Ridge dataset

Table 6. Comparison of average ERMSE values of different algorithms in the Jasper Ridge dataset

Table 7. Comparison of DSAD values of different algorithms in Urban dataset

Table 7. Comparison of DSAD values of different algorithms in Urban dataset

Table 8. Comparison of average ERMSE values of different algorithms in Urban dataset

Table 8. Comparison of average ERMSE values of different algorithms in Urban dataset

Table 1. Comparison of average D_SAD of various algorithms under the influence of Gaussian white noise

Table 1. Comparison of average D_SAD of various algorithms under the influence of Gaussian white noise

Table 2. E_RMSE comparison of various algorithms under the influence of Gaussian white noise

Table 2. E_RMSE comparison of various algorithms under the influence of Gaussian white noise

Table 3. Comparison of average D_SAD of various algorithms under the influence of salt and pepper noise

Table 3. Comparison of average D_SAD of various algorithms under the influence of salt and pepper noise

Table 4. E_RMSE comparison of various algorithms under the influence of salt and pepper noise

Table 4. E_RMSE comparison of various algorithms under the influence of salt and pepper noise

Table 5. Comparison of D_SAD values of different algorithms in the Jasper Ridge dataset

Table 5. Comparison of D_SAD values of different algorithms in the Jasper Ridge dataset

Table 6. Comparison of average E_RMSE values of different algorithms in the Jasper Ridge dataset

Table 6. Comparison of average E_RMSE values of different algorithms in the Jasper Ridge dataset

Table 7. Comparison of D_SAD values of different algorithms in Urban dataset

Table 7. Comparison of D_SAD values of different algorithms in Urban dataset

Table 8. Comparison of average E_RMSE values of different algorithms in Urban dataset

Table 8. Comparison of average E_RMSE values of different algorithms in Urban dataset