全变差稀疏约束深度非负矩阵分解高光谱遥感影像解混方法

赵文君，翟 晗，张洪艳

(1.武汉大学 测绘遥感信息工程国家重点实验室，湖北 武汉 430079；2.中国地质大学(武汉) 地理与信息工程学院，湖北 武汉 430074)

高光谱遥感结合了空间成像技术与光谱技术，能以纳米级的光谱分辨率获取地物近乎连续的光谱曲线，捕获地物间的细微差异，为地物精细化分类提供可能[1]，促进了遥感行业领域的进步与发展[2]，在农业、林业、食品安全、矿物勘探、环境监测等诸多领域得到广泛应用[3-6]。

由于成像光谱仪各波段获得的辐射能量有限，高光谱影像的空间分辨率通常较低[7]，再加上地物的复杂性，高光谱影像中存在大量混合像元，限制了高光谱遥感的行业应用水平。对此，高光谱解混从亚像元的角度出发，通过光谱混合模型将影像分解成纯净物质组分(端元)及其所占比例(丰度)来提高信息提取精度。

非负矩阵分解模型理论严谨，且易于与物理先验相结合，在高光谱解混中得到广泛应用。但由于目标函数的非凸性，导致解空间不唯一，降低了解混精度。故研究者们在目标函数中加入约束来提高算法精度。文献[8]将分段平滑度和稀疏性引入非负矩阵分解，同时保留了光谱数据的不连续性。文献[9]考虑本征流形结构，保持了原始影像和丰度图之间的紧密联系。文献[10]受到深度学习的启发，考虑了隐含信息的层次特征，解决了神经网络方法梯度扩散的问题。文献[11]引入自适应图来规范多层非负矩阵分解模型。

由于地物分布的复杂性，高光谱传感器接收到的信号通常是高度混合的。传统深度非负矩阵分解忽略了地物光谱的实际混合过程。为此，本文从地物光谱的传播出发，反向模拟光谱混合过程，提出面向端元矩阵的深度非负矩阵分解模型，并考虑端元的稀疏分布和空间分段平滑性，对丰度矩阵施加L1/2约束和全变差约束，提高了解混精度。

1 高光谱影像解混

1.1 线性光谱混合模型

线性光谱混合模型假设观测到的像元光谱是各端元光谱按一定比例线性混合而成的，如图1所示。线性混合模型可以用下式表示

X=AS+N

(1)

式中，X∈RB×N代表高光谱数据；B代表波段数；N代表像元数；A∈RB×P是端元光谱矩阵；P代表端元数；S∈RP×N是丰度矩阵；N∈RB×N代表噪声。

图1 线性光谱混合模型Figure 1. Explanation of linear spectral mixing model

由于地表地物分布复杂，到达高光谱传感器的信号通常是电磁波多次反射的结果，故高光谱传感器接收到的信号通常是邻近地物反射信号的深度融合。以电磁波信号在到达传感器前进行L次线性混合为例，深度线性光谱混合模型可用图2简要表示。

图2 深度线性光谱混合模型Figure 2. Explanation of deep linear spectral mixing model

则电磁波经过多次反射深度线性混合后得到的数据矩阵X如式(2)所示。

X=A1S1S2…SL+N

(2)

1.2 非负矩阵分解模型

非负矩阵分解[12](Nonnegative Matrix Factorization，NMF)将数据矩阵近似为两个非负矩阵的乘积，其代价函数如下式所示

(3)

式中，‖·‖F代表Frobenius范数。

1.3 深度非负矩阵分解模型

为研究高光谱遥感影像中的深度空谱信息，现有的深度非负矩阵分解从数学理论的角度出发，对高光谱数据进行深度分析，通过对分解过程中的系数矩阵(即丰度矩阵)进行深度分解，来探索高光谱遥感影像中隐含信息的特征。然而，若从数学理论层面对高光谱遥感影像解混进行分析，则会忽略高光谱遥感影像的成像过程及其所表征的物理意义。

从光谱混合的物理过程及其所表征的物理意义等方面考虑，由于地物光谱混合后的结果所表征的物理意义依旧为光谱，故对丰度矩阵进行深度分解是不合理的。这是由于其所表示的物理意义并非光谱，而是端元所占各像元的比例值。因此，本文从地物光谱的物理传播过程出发，对光谱深度线性混合的物理过程进行反向模拟，提出了面向端元矩阵的深度非负矩阵分解模型。

2 面向端元矩阵的深度非负矩阵分解

2.1 面向端元矩阵的深度非负矩阵分解模型

由于达到高光谱传感器的信号是电磁波在地表多次反射作用后的结果，故对原始数据进行一次分解所得到的端元矩阵不够纯净。针对该问题，本文从地物光谱的传播过程出发，对光谱深度线性混合过程进行反向模拟，提出了一种面向端元矩阵的深度非负矩阵分解模型。通过对每一层分解所得的端元矩阵进行分解，得到下一层的端元矩阵和丰度矩阵，直至得到最终的端元矩阵和丰度矩阵。

假定到达高空传感器的光谱从地物反射出发共经历了L次线性混合，设深度非负矩阵分解的层数为L，则高光谱遥感影像数据矩阵X的深度非负矩阵分解模型可表示为如图3所示的形式。

图3 深度非负矩阵分解模型示意图Figure 3. Explanation of deep NMF model

在第1层中，高光谱遥感影像数据矩阵X被分解为较纯净的端元矩阵A1及其丰度矩阵S1。在第l层中，上一层的端元矩阵Al-1被分解为更纯净的端元矩阵Al及其丰度矩阵Sl，则最终对高光谱遥感影像数据进行解混后所得到的纯净端元矩阵为A=AL，其丰度矩阵为S=SL…S2S1，则高光谱数据矩阵X可用下式表示

X=ALSL…S2S1+N

(4)

基于欧氏距离的代价函数可用式(5)进行表示。

(5)

本文提出的面向端元矩阵的深度非负矩阵分解模型包括预训练和微调两个阶段，如图4所示。预训练阶段通过最小化式(3)中的代价函数，用经典的非负矩阵分解对所有层进行逐层优化。微调阶段通过式(5)将预训练阶段所得的所有因素作为整体进行微调，从而减少总重建误差。

图4 深度非负矩阵分解流程图Figure 4. Illustration of deep NMF

2.2 基于全变差稀疏约束的深度非负矩阵分解模型

由于非负矩阵分解的代价函数非凸，存在大量极小值，因此有必要根据高光谱数据特点，对非负矩阵分解模型施加正则化约束。

研究表明，高光谱遥感影像的大多数像元仅由场景中某几个端元混合而成，其包含的端元数目通常远低于光谱库的端元总数，故丰度矩阵具有稀疏性。文献[13]证明了L1/2稀疏约束优于其它Lq正则化约束。基于此，本文将L1/2稀疏约束作为丰度矩阵的稀疏先验约束引入模型，得到预训练阶段和微调阶段的代价函数分别如式(6)和式(7)所示。

(6)

(7)

式中，λ为丰度稀疏约束项的系数，用于调节丰度稀疏约束的程度；‖·‖1/2代表L1/2稀疏约束，其计算式为

(8)

式中，si，j是影像中第i个端元占第j个像元的丰度；|·|代表绝对值。

高光谱影像具有空间平滑性。研究表明，全变差正则化能很好地提取影像的空间一阶邻域信息，可增加丰度图的分段平滑性并减少噪声的负面影响[14-16]。基于此，本文采用各向异性全变差正则化约束[17]来加强丰度图的分段平滑性，如下式所示。

(9)

综上，本文提出的全变差稀疏约束深度非负矩阵分解算法(Total Variation and Sparsity Regularized Deep Nonnegative Matrix Factorization，TVSDNMF)的预训练阶段和微调阶段的代价函数为

(10)

α‖Sl‖TV

(11)

式中，α为丰度全变差正则化项的系数，用于控制丰度图分段平滑程度。

3 实验结果与分析

3.1 实验设置

本文通过模拟实验和真实实验分别验证所提出的全变差稀疏约束的深度非负矩阵分解算法对高光谱遥感影像解混的有效性。本文选择顶点成分分析-全约束最小二乘(Vertex Component Analysis-Fully Constrained Least Squares，VCA-FCLS)[18]、L1/2稀疏约束非负矩阵分解(L1/2Sparsity-constrained Nonnegative Matrix Factorization，L1/2-NMF)[13]、图正则化L1/2稀疏约束非负矩阵分解(Graph Regularization of Nonnegative Matrix Factorization，GNMF)[9]、多层非负矩阵分解(Multilayer Nonnegative Matrix Factorization， MLNMF)[19]和全变差稀疏约束的深度非负矩阵分解(Sparsity-constrained Deep Nonnegative Matrix Factorization with Total Variation， SDNMF-TV)[10]5种高光谱解混算法作为对比，并采用光谱角距离和均方根误差分别对光谱提取及丰度估计结果进行定量评价。

(12)

(13)

3.2 模拟实验

本文采用的模拟数据集由美国地质调查局光谱库中的数据生成。该库包含近500种典型矿物，有224个波段，均匀分布在0.4～2.5 μm之间。本文选择光卤石、胺基明矾石、黑云母和阳起石的188个低噪声波段进行模拟实验，其光谱曲线如图5所示。模拟数据集由48×48像元组成[20]。数据是使用线性光谱混合模型生成的，同时在各模拟像元中施加丰度和唯一约束。在模拟数据中，存在纯净区域和混合区域两种不同的区域，第1行为纯净区域，仅由一种端元生成，第2～4行为混合区域，分别使用两个端元、3个端元和4个端元混合生成。各区域以方形分布，其具体分布如图6所示。

图5 模拟实验端元光谱图Figure 5. Spectral signatures of synthetic experiments

图6 模拟实验丰度图Figure 6. Abundance maps of synthetic experiments

表1给出了使用不同解混方法得到的模拟数据集平均光谱角距离值和均方根误差的定量评价，其中加粗字体为每列最优值。从表中可以看到，TVSDNMF的平均光谱角距离值远低于其他方法。相比于面向丰度矩阵的深度非负矩阵分解方法，平均光谱角距离减小了0.007 8，均方根误差降低了0.006 8，证明了面向端元矩阵的深度非负矩阵分解方法的优越性与合理性。

表1 不同算法模拟实验定量评价结果

图7示出了使用不同解混方法提取到的端元光谱，其中实线为端元真实光谱，虚线为端元估计光谱。可以看到，除了MLNMF，其余各方法端元提取结果精度较高。相比于面向丰度矩阵的深度非负矩阵分解方法，面向端元矩阵的深度非负矩阵分解方法的吻合度更高，再次证明了所提出的面向端元矩阵的深度非负矩阵分解方法的优越性。

(a)

图8则显示了使用不同解混方法提取到的端元1所对应的丰度图。可以看出，两种深度非负矩阵分解的丰度图的噪声点相较其他算法更少。通过对两种方法进行比较，可以发现，所提出的TVSDNMF算法得到的丰度图的背景噪声更小。

(a) (b) (c)

3.3 真实实验

本文采用了Samson高光谱数据集进行真实实验，进一步验证了所提算法的有效性。Samson数据集包含了952×952个像元，每个像元有156个波段，波段范围为0.401～0.889 μm，其光谱分辨率高达3.13 nm。由于原始影像范围较大，这会使得计算时间过长，故本文参考之前的研究选择了(252，332)处95×95个像元的区域来减少计算量，因为此区域的数据不会由于空白光谱波段或严重噪声的影响而降低精度[21]。该选区影像中主要含有3个端元，分别是土壤、树木和水，所用选区影像图如图9所示。

图9 Samson高光谱数据集Figure 9. Samson dataset used in real data experiment

表2给出了使用不同解混方法得到的Samson数据集平均光谱角距离值和均方根误差的定量评价， 其中加粗字体为每列最优值。从表中可以看到，TVSDNMF方法的平均光谱角距离值远低于其他方法，相比于面向丰度矩阵的深度非负矩阵分解方法，其平均光谱角距离即均方根误差均有所降低，这证明了所提方法的优越性。

图10示出了使用不同解混方法提取的端元光谱曲线，其中实线为端元真实光谱，虚线为端元估计光谱。从图中可以看到，两种深度非负矩阵分解方法提取到的端元估计光谱与端元真实光谱更为一致。相较于SDNMF-TV，TVSDNMF估计所得的水的光谱曲线与真实光谱更为贴合。

(a)

图11～图13显示了使用不同解混方法提取到的不同端元所对应的丰度图。可以看出，整体而言，两种深度非负矩阵分解方法SDNMF-TV与TVSDNMF估计所得的丰度图与真实值更为接近。在土壤的丰度图中，除了两种深度非负矩阵分解方法，其余方法所得丰度值均较小。相较两种深度非负矩阵分解方法，本文所提出的TVSDNMF丰度估计结果更为细致，证明了所提方法的优越性。

表2 Samson数据集不同算法的定量评价结果

(a) (b) (c) (d) (e) (f) (g) 图11 不同方法提取的Samson数据集端元“水”的丰度图(a)地面真实值 (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMFFigure 11. Abundance map of water in Samson dataset extracted by different methods(a)Ground truth (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMF

(a) (b) (c) (d) (e) (f) (g) 图12 不同方法提取的Samson数据集端元“树”的丰度图(a)地面真实值 (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMFFigure 12. Abundance map of tree of Samson dataset extracted by different methods(a)Ground truth (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMF

图13 不同方法提取的Samson数据集端元“土壤”丰度图(a)地面真实值 (b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e) MLNMF (f)SDNMF-TV (g)TVSDNMFFigure.13 Abundance map of soil of Samson dataset extracted by different methods(a)Ground Truth(b)VCA-FCLS (c)L1/2-NMF (d)GNMF (e)MLNMF (f)SDNMF-TV (g)TVSDNMF

4 结束语

本文针对现有深度非负矩阵分解方法未考虑地物光谱实际混合过程的问题，从地物反射电磁波的传播过程出发，对光谱深度线性混合过程进行反向建模，提出了面向端元矩阵的深度非负矩阵分解模型。此外，本文引入L1/2稀疏约束和全变差正则化来增强丰度稀疏性和空间平滑性，进一步提高了解混精度，建立了全变差稀疏约束的深度非负矩阵分解模型。与5种解混方法(VCA-FCLS、L1/2-NMF、GNMF、MLNMF和SDNMF-TV)的对比结果表明，本文所提出的全变差稀疏约束的深度非负矩阵分解方法能更好地提取端元，具有良好的性能。

今后的工作将集中于将本文模型扩展到能够应对混合噪声(稀疏噪声和高斯噪声)的情形，以提高模型对噪声的鲁棒性。