高光谱遥感利用数百个波段采集地表地物从可见光至近红外的光谱响应信息，能够利用分类区分具有细微光谱特征差异的地物^[1-3]，分类结果对植被覆盖制图和海岸带监测等应用影响重大^[4-5]。但是，高光谱波段众多且相关性强，导致分类时存在明显的信息冗余且计算量较大^[6]。此外，高光谱分类存在维数灾难问题，需要巨大的训练样本来得到高精度的结果^[7]。因此，可以采用波段选择来选取波段子集，降低原始数据计算量并得到较高精度^{[3, 8]}。

国内外学者针对高光谱波段选择提出了许多方法，包括基于排序策略的最大方差主成分分析（maximum variance principal component analysis, MVPCA）^[9]与限制能量最小的限制波段选择方法^[10]，基于聚类策略的层次聚类法^[11]与仿射传播法^[12]，以及基于搜索策略的线性预测搜索法^[13]与粒子群优化搜索法^[14]等。其后，矩阵计算策略逐渐被用于解决高光谱波段选择问题，充分考虑高光谱数据的潜在结构（如稀疏、低秩和非负等），将波段选择转换为矩阵分解的目标优化问题来求解。目前有学者在该方面进行了研究，如稀疏非负矩阵分解（sparse nonnegative matrix factorization, SNMF）将高光谱波段矩阵分解为字典和稀疏系数矩阵，利用聚类稀疏系数的列向量来选择最佳波段^[15]；稀疏支持向量机方法基于二值分类器得到各个波段对应的权重系数，利用权重系数差异来选取重要波段^[16]；改进的稀疏子空间聚类方法通过整合稀疏表达和子空间聚类方法，利用稀疏系数向量来构建相似性矩阵，采用谱聚类方法来得到合适波段^[6]；多任务稀疏追踪将波段选择转换为多任务稀疏追踪学习问题，采用进化策略来选取波段^[17]；差异加权自表达考虑波段的差异信息来构建差异性正则项，改进自表达模型来选取较大差异的波段^[8]；低阶表达方法利用排序波段的表达系数来选取合适波段以满足分类需求^[18]。

不同于上述方法，本文提出了可分离非负矩阵分解(separable nonnegative matrix factorization, SepNMF）方法来研究高光谱影像的波段选择问题。利用美国印第安纳派恩斯市和某城市共两个数据集进行实验分析，验证提出的SepNMF方法的有效性。

1 可分离非负矩阵分解实现波段选择 1.1 高光谱影像的可分离非负矩阵分解模型

假设高光谱波段矩阵为$\mathit{\boldsymbol{Y}} = [{\mathit{\boldsymbol{y}}_i}]_{i = 1}^N \in {{\bf{R}}^{D \times N}}$，其中D为像素个数，N为波段数，y_i为任一波段对应的向量。非负矩阵分解模型为：

$\mathit{\boldsymbol{Y}} = \mathit{\boldsymbol{WH}} $

(1)

该模型将高光谱波段矩阵Y分解为非负字典矩阵W∈R^D×r和非负系数矩阵H∈R^r×N的乘积，以逼近原始的影像数据^[19-20]，其中r为W和H的秩，且r≤N < D。非负矩阵分解具有许多特性，如聚类、稀疏和低秩等，然而W和H未知，导致上述问题非凸而较难得到全局最优解。

可分离特性能够解决上述病态问题^[21]。可分离特性认为^[22]：假如存在一个基数为r的高维数据Y的索引子集k满足Y=Y(:, k_i)H，那么Y是r分离的。其中，“：”表示选中矩阵中某列全部元素，Y(:, k_i)为水平排列的索引子集k的原始波段矩阵Y的子集，$\forall {k_i} \in \left[{1, 2 \cdots N} \right]$且|k|₀=r。

将可分离特性应用到Y中，假设W=Y(：，k_i)，高光谱的可分离非负矩阵分解模型可表达为：

$\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Y}} = \mathit{\boldsymbol{WH}}, \forall {k_i} \in \left[ {1, 2 \cdots N} \right], }\\ {\quad {\rm{}}{{\left| k \right|}_0} = r{\rm{}}{\rm{}}r \ll {\rm{min}}\left( {D, N} \right)} \end{array} $

(2)

式中，k为波段的索引子集且数量为r。可分离特性能够保证高光谱波段矩阵Y估计得到唯一的波段子集Y(, :, k_i)，且剩余的波段能够被Y(:, k_i)进行非负线性表达。可分离矩阵分解求解解释为寻找包围其对应高维数据点的凸圆锥体的极端射线^[23]。

1.2 可分离非负矩阵分解模型的求解

式（2）中，假如将高光谱波段矩阵Y的每一列进行归一化处理，则可以得到：

$\mathit{\boldsymbol{Y}} = \mathit{\boldsymbol{WH}} \equiv \mathit{\boldsymbol{YD}}_\mathit{\boldsymbol{Y}}^{ - 1} = \mathit{\boldsymbol{WD}}_\mathit{\boldsymbol{W}}^{ - 1}\left( {{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{W}}}\mathit{\boldsymbol{HD}}_\mathit{\boldsymbol{W}}^{ - 1}} \right) $

(3)

式中，

可以看出，D_Y^-1和WD_W^-1每一列元素和为1。同时，由可分离特性可知，高光谱的波段矩阵可以通过优化排列，使得前r列分别对应于W的列，即Y=WH=W[I_rH′]^T，其中I_r为大小为r的单位矩阵，H′∈R^(N-r)×r为系数矩阵，其每一列所有元素和小于等于1。考虑到采集误差，假设Y′=Y+E，其中E为误差项，Y满足可分离特性且每一列元素和为1。

文献[24]证明，如果索引i使得Y′(：，i)满足条件argmax‖Y′(：，i)‖₂，则存在j∈k使得W(：，j)满足：

$\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Y}}{\rm{'}}\left( {:, i} \right) - \mathit{\boldsymbol{W}}{{\left( {:, j} \right)}_2} \le O\left( {\frac{L}{\mu }\frac{\varepsilon }{{\sigma _r^2\left( \mathit{\boldsymbol{W}} \right)}}} \right), }\\ {1 \le i \le N, {\rm{}}1 \le j \le r} \end{array} $

(4)

式中，μ和L分别为函数f=‖·‖₂的凸参数和利普希茨连续梯度常数；‖·‖₂的全局最小值为f(0)=0；ε为误差项的二阶范数阈值，即‖E(：，i)‖≤ε；σ_r(W)为矩阵W的第r个奇异值，O(·)表示和括号里的参数成比例的一个数值。本文采用迭代优化策略，假定第t次迭代中Y′的残差矩阵R^(t)= [W′^(t)Q]H+E^(t)，其中Q∈R^D×(r-t)，可分离非负矩阵分解模型式（2）可转换为优化问题式（5）：

${i^{\rm{*}}} = {\rm{argmax}}\left\| {{\mathit{\boldsymbol{R}}^{\left( t \right)}}{{\left( {:, i} \right)}_2}} \right\|, 0 \le t \le r $

(5)

可以看出，在第t次迭代中，argmax‖R^(t)(：，i)‖₂能够保证找到j∈κ来满足式(4)的条件。本文采用迭代投影算法^[18]来求解目标优化问题，即式（5）。该方法利用R⁽⁰⁾=Y′来初始化残差矩阵R⁽⁰⁾，采用W(：，t)=Y′(：，i^*)来更新波段子集，并通过${\mathit{\boldsymbol{R}}^{\left({t + 1} \right)}} = \left({\mathit{\boldsymbol{I}} - \frac{{\mathit{\boldsymbol{R}}\left({:, {i^{\rm{*}}}} \right)\mathit{\boldsymbol{R}}{{(:, {i^{\rm{*}}})}^{\rm{T}}}}}{{\mathit{\boldsymbol{R}}\left({:, {i^{\rm{*}}}} \right)_2^2}}} \right){\mathit{\boldsymbol{R}}^{\left(t \right)}}$来迭代更新下一步的残差矩阵。在迭代次数t=r后，可估计得到合适的波段子集W=Y′(：，k)，从而进一步优化求解目标函数式(6)来估计得到系数矩阵Ĥ：

$\mathit{\boldsymbol{\hat H}} = \mathop {{\rm{argmin}}}\limits_{\mathit{\boldsymbol{H}} > 0} \left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{WH}}} \right\|_\mathit{\boldsymbol{F}}^2 $

(6)

1.3 SepNMF实现波段选择的流程

利用可分离非负矩阵分解方法SepNMF实现高光谱波段选择的流程如图 1所示，包含以下步骤：

1）将三维高光谱影像立方体转换为二维矩阵Y′，其中D为Y′的矩阵列数，N为Y′的矩阵行数。

2）利用式(2)构建高光谱影像的可分离非负稀疏矩阵分解模型，并将式(2)转换为寻找最大化残差矩阵列的目标优化问题，即式(5)。

3）利用迭代投影算法依次更新误差矩阵R^(t+1)和W(：，t)。在迭代次数t=r时，估计得到合适的波段子集W=Y′(：，k)。

4）求解目标函数式(6)得到系数矩阵的最优值Ĥ。

2 波段选择实验与效果分析 2.1 实验数据获取与预处理

美国印第安纳派恩斯市影像来自美国普渡大学遥感应用实验室。影像由美国喷气推进实验室成像光谱仪于1992年6月12日采集得到，空间分辨率为20 m，光谱分辨率为10 nm，光谱区间为200~2 400 nm，预处理后波段数为200个。图 2为覆盖美国印第安纳州西拉法叶地区西部8 km（派恩斯市）的一小块区域，包含145×145像素。图 2中共包含16类主要地物，其地物分布的样本信息见表 1。

美国某城市高光谱影像为从美国陆军地理空间中心获取的HYDICE影像。影像采集于1995年10月，空间分辨率为2 m，光谱分辨率为10 nm。影像大小为307×307像素，覆盖美国德克萨斯州科帕拉斯区域（靠近胡德堡），如图 3所示。对原始的210个波段数据进行预处理，移除低噪比波段，剩余162个波段，包含22类主要地物，各地物的真实分布样本信息如表 2所示。

2.2 SepNMF波段选择效果分析

利用公开高光谱影像（印第安纳派恩斯市和某城市）设计实验来验证SepNMF波段选择的效果。实验采用支持向量机（support vector machine, SVM）分类器和k-近邻（k-nearest neighbor, kNN）分类器来分类波段子集，利用总体分类精度（overall classification accuracy, OCA）来评价分类效果。SVM采用径向基核函数，其方差和惩罚因子通过交叉验证获得，kNN分类器的邻域大小为1。对比方法包括MVPCA ^[9]、基于稀疏的波段选择法（sparse based band selection, SpaBS）^[25]、SNMF^[15]和基于快速密度峰值聚类法（fast density-peak based clustering, FDPC）^[26]。

2.2.1 不同波段选择方法的定量评价对比

本实验对比不同波段选择方法得到的波段子集的定量评价结果。参照文献[6]，实验采用平均信息熵（average information entropy, AIE）来定量评价波段子集的信息量及光谱信息的丰富程度，采用平均相关系数（average correlation coefficient, ACC）来定量评价各波段的相关关系，采用平均相对熵（average relative entropy, ARE）来定量评价各个波段的差异及用于分类的可分性。印第安纳派恩斯市和某城市数据中，k的大小分别选取为30和40，SpaBS的迭代次数设置为5次。印第安纳派恩斯市数据中，SNMF的正则化因子分别设定为3.0和0.1；某城市数据中，SNMF方法中正则化因子分别设定为3.5和0.05。

表 3对比了不同波段方法的定量评价结果。SepNMF波段子集的AIE最高，ACC最低，说明波段信息量最大、相关性最小。SepNMF的ARE表现较好，在印第安纳派恩斯市数据中略低于FDPC，但优于其他4种方法。SNMF和FDPC方法的定量评价结果总体较为接近，MVPCA和SpaBS方法的定量评价结果低于其他3种方法，在5种方法中表现欠佳。

2.2.2 不同波段数量下的分类结果对比

对比分析不同波段子集大小条件下SepNMF和其他几种方法在印第安纳派恩斯市和某城市数据集的分类性能。实验中，印第安纳派恩斯市影像中波段数的选择区间为2~44，步长为2；某城市影像中波段数的选择区间为2~50，步长为2。实验选定每一类地物中20%的真实样本作为训练样本来进行分类实验，其余地物样本作为测试样本。

图 4为在不同子集大小条件下，各种波段选择方法在印第安纳派恩斯市和某城市影像集上得到的OCA曲线图。可以看出，随着波段数量的增加，所有方法的OCA整体上升。SepNMF方法得到的OCA曲线结果最优，分类精度整体优于其他几种方法。相比而言，SpaBS的OCA曲线整体低于MVPCA和其他几种方法（SNMF、FDPC和SepNMF）。FDPC的OCA曲线与SNMF较为接近或略低，但优于MVPCA和SpaBS方法。图 5和图 6展示了以上方法在选定波段子集数量条件下的SVM分类结果。其中，印第安纳派恩斯市数据和某城市数据中，波段子集k的大小分别为30和40。可以看出，SepNMF的分类图结果优于其他方法，进一步验证了图 4的结论。因此，在不同波段数量条件下，SepNMF的OCA最高，优于其他几种波段选择方法。

2.2.3 不同训练样本下的分类结果对比

此实验对比分析不同训练样本大小条件下SepNMF和其他4种对比方法（MVPCA、SpaBS、SNMF和FDPC）的分类结果。其中，印第安纳派恩斯市和某城市影像的训练样本随机采样于地面地物样本集，采样比率的选择集合为为{0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5}。印第安纳派恩斯市影像中，波段数量k设定为30；某城市影像中，波段数量k设定为40。

图 7为在不同地物训练样本大小条件下，SepNMF和其他几种方法在印第安纳派恩斯市和某城市影像中得到的OCA曲线图。可以看出，随着单一地物训练样本数量的增加，不同波段选择方法的OCA明显提升。其中，SepNMF方法得到的OCA曲线最优，在印第安纳派恩斯市和某城市影像中均优于其他方法的OCA曲线。SNMF和FDPC的OCA曲线结果较为接近；当训练样本较大时，SNMF的OCA结果优于FDPC。MVPCA的OCA曲线低于SNMF、FDPC和SepNMF，但明显优于SpaBS。因此，SepNMF在不同训练样本条件下得到的OCA结果优于其他几种波段选择方法。

2.2.4 波段选择方法的计算时间对比

此实验对比SepNMF和其他4种方法在不同波段数量条件下的计算效率。表 4为不同波段数量条件下，5种波段选择方法在印第安纳派恩斯市和某城市影像集上得到的计算时间。实验环境为戴尔Win10台式计算机，主要配置为Inter Core (TM) i7-6700 CPU 3.40 GHz和32 GB内存，所有对比方法的代码执行环境为Matlab 2016a。可以看出，随着波段数量的增加，各种方法的计算时间呈上升趋势。与SNMF和SpaBS相比，FDPC、SepNMF和MVPCA的计算时间随波段数量的增加上升幅度较小。这5种方法中，SpaBS计算效率最低，计算时间最长。SNMF的计算速度优于SpaBS，但低于其他方法。FDPC的计算速度低于SepNMF和MVPCA，MVPCA的计算速度最快。因此，5种方法的计算效率排名为：MVPCA > SepNMF > FDPC > SNMF > SpaBS。

3 结语

本文提出利用SepNMF方法来解决高光谱波段选择问题，假设波段集合具有可分离结构，将波段选择转换为寻找非负线性表达其他波段的代表列，采用迭代投影方法来选取适合的波段子集。实验表明，SepNMF方法在不同波段数量和不同训练样本条件下的OCA优于其他4种方法，能够得到更好的分类结果。同时计算效率较好，在5种波段选择方法中排列第2。然而，本文没有考虑波段间的细微差异可能对结果产生的影响，在后续的研究中，将考虑波段差异信息权重矩阵，细致分析噪声和异常值对波段选择的影响。