依赖于成像光谱仪，高光谱遥感技术能够以数十至数百个窄波段来收集地物可见光至近红外的光谱响应信息(童庆禧 等，2016；张良培和李家艺，2016)。在此基础上，通过分类来区分具有细微光谱差异的地物，对后续的植被覆盖制图、海岸带监测、矿物制图和军事目标识别等应用关系重大(杜培军 等，2016)。然而，高光谱影像数据由于波段众多，导致分类存在明显的“维数灾难”问题，需要巨大的训练样本来保证高精度的分类结果(Pal和Foody，2010)。另一方面，波段众多且相关性强，导致高光谱影像内部信息冗余较高，数据处理的计算量较大。因此，通常情况下，可以采用波段选择来进行高光谱影像数据降维，通过选取波段子集来降低数据计算量并保留重要光谱信息，然后再进行后续的高光谱影像分类处理(Sun 等，2016a, 2015b)。

近年来，稀疏表达理论的兴起，为高光谱影像的波段选择提供了新的研究思路。稀疏表达理论认为任一正常信号可以在特定字典的条件下表达为一个稀疏向量，其中少数非零元素能够很好继承原始信号的特性(Donoho，2006)。稀疏表达理论目前已经广泛应用于高光谱影像的混合像元分解(Du 等，2016)、分类(Du 等，2015)和异常探测(Zhang 等，2016)等领域。高光谱影像的稀疏系数向量能够凸显原始影像内部的一些潜在结构信息如聚类和潜在子空间等，便于实现影像的波段选择(Sun 等，2015a；张良培和李家艺，2016)。Chepushtanova等人(2014)提出稀疏支持向量机方法，利用二值支分类器来训练得到高光谱影像的各波段对应的权重，通过零和非零权重的差异来选取重要波段并得到波段子集。施蓓琦等人(2013)分析稀疏系数矩阵的聚类结构特征，利用欧氏距离度量和相对熵距离度量约束将高光谱数据矩阵分解为低阶字典矩阵和稀疏系数矩阵，根据稀疏系数矩阵中对应的系数权重大小确定该波段的聚类，最终实现波段的有效选取。Sun等人(2015b)提出改进稀疏子空间聚类方法，整合稀疏表达和子空间聚类方法，通过稀疏系数向量来构建相似性矩阵并采用谱聚类来得到合适的波段子集。最近，Sun等人(2016a)整合差异性正则项和稀疏自表达模型来构建加权差异性稀疏自表达的波段选择方法，通过非线性求解稀疏表达系数来实现波段子集的优化选择。

不同于以上方法，本文提出加权概率原型分析WPAA (Weighted Probabilistic Archetypal Analysis)方法来研究高光谱影像的波段选择问题。通过两个高光谱影像公开数据集包括Urban和Pavia University(PaviaU)进行实验分析，验证本文提出WPAA方法的有效性。

原型分析理论由Culter和Breiman(1994)提出，目前广泛应用于寻找高维数据集的代表性样本点。对于高光谱影像而言，该理论认为数据集中任一波段向量能够由其目标波段子集进行稀疏表达的线性组合来重构得到；同时，任一代表性波段可以视为由所有原始波段向量进行稀疏线性组合得到，即波段子集对应的数据点位于原始波段对应的高维数据集合所构成的凸包的边界上(Sun 等，2016a)。因此，对于高光谱影像，原型分析理论认为选取波段子集为寻找原始所有波段向量对应的高维数据点构成的最小凸包的原型点(Mørup和Hansen，2012)。假设高光谱影像构成波段集合为 $\mathit{\boldsymbol{Y}} = \left\{ {{\mathit{\boldsymbol{y}}_j}} \right\}_{j = 1}^N \in {{\bf{R}}^{D \times N}}$ ，其中D为像素个数，N为波段数量， $\mathit{\boldsymbol{Z}} = \left\{ {{{\textit{z}}_l}} \right\}_{l = 1}^k \in {{\bf{R}}^{D \times k}}$ 为待求的波段子集，其中k为波段子集的大小。高光谱影像的原型分析模型如下(Sun 等，2016b)：

式中， $\mathit{\boldsymbol{A}}{\rm{ = }}\left\{ {{\mathit{\boldsymbol{a}}_i}} \right\}_{i = 1}^N$ 和 $\mathit{\boldsymbol{B}} = \left\{ {{\mathit{\boldsymbol{b}}_j}} \right\}_{j = 1}^k$ 为稀疏系数矩阵，E为由于模型误差和随机噪声导致的误差项。 ${\mathit{\boldsymbol{a}}_i} \geqslant 0$ 和 $\left| {{\mathit{\boldsymbol{a}}_i}} \right| = 1$ 限制是保证波段子集用以重构原始波段集合中任一波段的概率为1且任一波段的概率权重都大于0。 ${\mathit{\boldsymbol{b}}_{jj}} = 0$ ， ${\mathit{\boldsymbol{b}}_j} \geqslant 0$ 和 $\;\left| {{\mathit{\boldsymbol{b}}_j}} \right| = 1$ 是保证任一原始波段用以重构波段子集的概率权重都大于0且重构概率为1。

常规的原型分析模型通过选取原型点来估计最佳波段子集时，并未考虑到高光谱影像不同波段的差异性导致的原型点选取时的优先度问题。波段选择的目的为选择信息量大且相关性较小的波段，因此本文采用综合差异度量CDM(Composite Dissimilarity Measure)来描述不同波段之间的综合差异，改进原型分析模型用于高光谱影像的波段选择结果(Sun 等，2016a)。综合差异性度量整合光谱信息差异性和波段相关性指标来综合评价不同波段的差异大小。假设任意两个波段向量 ${\mathit{\boldsymbol{y}}_i} = $ $ {({a_{i1}}, \cdot \cdot \cdot {a_{il}}, \cdot \cdot \cdot ,{a_{iD}})^{\rm{T}}} \in {{\bf{R}}^{D \times 1}}$ 和 ${\mathit{\boldsymbol{y}}_j} = {({a_{j1}}, \cdot \cdot \cdot, {a_{jl}}, \cdot \cdot \cdot ,{a_{jD}})^{\rm{T}}} \in $ $ {{\bf{R}}^{D \times 1}}$ ，其中 ${a_{il}}$ 和 ${a_{jl}}$ 为像素 ${y^l}$ 在第i和第j个波段影像上对应的光谱值，综合差异性度量CDM的定义如下(Sun 等，2016a)

式中， ${\rm{KLD}}(\left. {{\mathit{\boldsymbol{y}}_i}} \right\|{\mathit{\boldsymbol{y}}_j}) = \sum\limits_{l = 1}^N {{{\hat a}_{il}}\lg ({{\hat a}_{il}}/{{\hat a}_{jl}})} $ 是波段 ${\mathit{\boldsymbol{y}}_i}$ 对于 ${\mathit{\boldsymbol{y}}_j}$ 的相对信息熵， ${\hat a_{il}}$ 为正则化的光谱响应值， ${\bar y_i}$ 和 ${\bar y_j}$ 为波段 ${\mathit{\boldsymbol{y}}_i}$ 和 ${\mathit{\boldsymbol{y}}_j}$ 的平均光谱值。在此基础上，可以得到波段差异性矩阵 $\mathit{\boldsymbol{W}} = \left\{ {{\bf{CD}}{{\bf{M}}_{ij}}} \right\} \in {{\bf{R}}^{N \times N}}$ 来描述波段之间的综合差异性。通过在 $\mathit{\boldsymbol{Z}} = \mathit{\boldsymbol{YB}}$ 中附加差异性矩阵作为权重系数来添加波段的差异性信息至稀疏系数矩阵B，得到加权原型分析模型公式如下

式中，稀疏系数矩阵B和A均未知，是典型的盲信号矩阵分解问题。差异性矩阵通过权重稀疏来综合考虑不同波段的信息量和相关性，避免加权原型分析选择较为相近的波段，更有利于选择信息量大且相关性较小的波段子集。贝叶斯框架考虑光谱的空间变异性，能够将矩阵分解问题(5)转换为基于概率模型的贝叶斯预测问题。贝叶斯框架较为灵活，能够整合波段和光谱信号的先验概率分布知识，同时附加各类范数约束条件如稀疏系数矩阵B和A的非负和一阶范数为1的限制条件来保证优化结果能够得到合理解释(Bioucas-Dias 等，2012；He 等，2016)。得到加权概率原型分析模型

式中， $\mathit{\boldsymbol{\varTheta }} = \mathit{\boldsymbol{YW}}$ ，P(·)是概率密度函数， $P(\mathit{\boldsymbol{\varTheta }})$ 和 $P(\mathit{\boldsymbol{Y}})$ 为常数。式(6)的联合概率密度最大化问题可以转换为最大后验估计问题

式中，lg(·)为对数运算符， $P(\mathit{\boldsymbol{Y}}\left| {\mathit{\boldsymbol{\varTheta }},\mathit{\boldsymbol{B}},\mathit{\boldsymbol{A}}} \right.)$ 依赖于高光谱成像的观测模型的假设条件。

考虑到高光谱遥感的光学成像的量子特性，本文采用泊松分布来描述光谱的概率密度分布 $P(\mathit{\boldsymbol{Y}}\left| {\mathit{\boldsymbol{\varTheta }},\mathit{\boldsymbol{B}},\mathit{\boldsymbol{A}}} \right.)$ (Patil 等，2016)。考虑到稀疏系数矩阵B和A的非负和一阶范数为1中各元素非负且一阶范数为1，因此采用狄利克雷函数来描述概率密度分布，即 ${\mathit{\boldsymbol{a}}_i} \sim Dir(0,1)$ 和 ${\mathit{\boldsymbol{b}}_j} \sim Dir(0,1)$ (Dobigeon 等，2014)。进一步，最大后验估计问题式(7)转换为

假设 $ < \cdot > $ 代表矩阵的内积， $ \circ $ 代表矩阵的阿达玛积，上述目标优化问题式(8)可以转换为

式中， $C$ 代表 $\lg {(\mathit{\boldsymbol{Y}})_{ij}}!$ 产生的常数项，1为所有元素为1的矩阵。可以采用迭代求解的策略，在t次迭代过程中，固定稀疏系数矩阵B或A的情况下，式(9)可以分解为两个凸优化子问题

本文采用交替方向乘积方法(Sun 等，2015b)来分别求解上述两个凸优化子问题，通过引入适当的辅助变量来增强上述约束条件，依次最小化主要变量和拉格朗日乘积构成的目标函数方程来得到局部最优解。当目标残差满足 ${\left\| {{\mathit{\boldsymbol{B}}^{(t + 1)}} - {\mathit{\boldsymbol{B}}^{(t)}}} \right\|_\infty } \leqslant $ $ \varepsilon $ ， ${\left\| {{\mathit{\boldsymbol{A}}^{(t + 1)}} - {\mathit{\boldsymbol{A}}^{(t)}}} \right\|_\infty } \leqslant \varepsilon $ ， ${\left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{\varTheta }}{\mathit{\boldsymbol{B}}^{(t + 1)}}{\mathit{\boldsymbol{A}}^{(t + 1)}}} \right\|_\infty } \leqslant \varepsilon $ 或迭代次数t大于设定阈值时，停止迭代并得到最终的稀疏系数矩阵 $\hat{B} = {\mathit{\boldsymbol{B}}^{(t + 1)}}$ 和 $\hat{A} = {\mathit{\boldsymbol{A}}^{(t + 1)}}$ 。

基于加权概率原型分析WPAA的高光谱影像波段选择的流程如图1。

流程中包含以下步骤：

(1) 将3维的高光谱影像立方体转换为2维波段向量数据集Y，以2维矩阵形式表示，其中波段数D为Y对应的矩阵的列数，像素个数N为Y的行数；

(2) 利用式(6)构建加权概率原型分析模型，并通过联合概率密度最大化来得到基于最大后验估计函数的目标优化问题式(7)；

(3) 利用迭代优化的策略来分解问题式(7)为两个凸优化子问题式(10)和式(11)，采用交替方向乘积方法来求解得到局部最优的稀疏系数矩阵 $\hat{B} = {\mathit{\boldsymbol{B}}^{(t + 1)}}$ 和 $\hat{A} = {\mathit{\boldsymbol{A}}^{(t + 1)}}$ ；

(4)利用稀疏重构表达式 $\hat{Z} = {\mathit{\boldsymbol{YW}}}{\hat {B}}$ 得到待求的波段子集 $\hat{Z}$ 。

图 1 WPAA用于高光谱影像波段选择的技术流程 Figure 1 The procedure of band selection with WPAA on Hyperspectral imagery

PaviaU数据来自西班牙巴斯克大学计算智能课题组(http://www.ehu.es/ccwintco/index.php/Hyper spectral_Remote_Sensing_Scenes [2016-10-11])，影像覆盖帕维亚大学区域，共103波段，空间分辨率为1.3 m(图2)。影像为较大数据集中的一部分，包含350×340像素，波段数为103，包含9类地物(包括阴影)。数据中各地物样本的真实空间分布信息如表1所示。Urban数据从美国陆军地理空间中心获取的HYDICE高光谱影像数据(http://www.erdc.usace.army.mil/Media/FactSheets/FactSheetArticleView/tabid/9254/Article/610433/hypercube.aspx[2016-10])。数据采集于1995年10月，空间分辨率为2 m，光谱分辨率为10 nm。影像大小为307×307像素，覆盖美国德克萨斯州科帕拉斯区域(靠近胡德堡)(图3)。对原始的210波段数据进行预处理，移除低噪比波段，剩余162波段，包含22种主要地物，各地物的真实分布样本信息如表2。

图 2 PaviaU数据 Figure 2 The image of PaviaU data

图 3 Urban数据 Figure 3 The image of Urban data

表 1(Table 1)

表 1 PaviaU数据的地物真实样本信息 Table 1 Ground truth of all classes of ground objects on PaviaU dataset 类号 类名 样本数 类号 类名 样本数 1 Asphalt 4195 6 Bare Soil 5029 2 Meadows 2185 7 Bitumen 1330 3 Gravel 2099 8 Self-Blocking Bricks 2347 4 Trees 1550 9 Shadows 929 5 Painted metal sheets 1345 总数 20289

表 1 PaviaU数据的地物真实样本信息 Table 1 Ground truth of all classes of ground objects on PaviaU dataset

表 2(Table 2)

表 2 Urban数据的地物真实样本信息 Table 2 The ground truth of samples in each class for Urban dataset 类号 类名 样本数 类号 类名 样本数 1 AsphaltDrk 85 13 Roof03LgtGray 35 2 AsphaltLgt 58 14 Roof04DrkBrn 84 3 Concrete01 124 15 Roof05AChurch 85 4 VegPasture 236 16 Roof06School 64 5 VegGrass 127 17 Roof07Bright 72 6 VegTrees01 263 18 Roof08BlueGrn 45 7 Soil01 113 19 TennisCrt 96 8 Soil02 53 20 ShadedVeg 40 9 Soil03Drk 59 21 ShadedPav 64 10 Roof01Wal 118 22 VegTrees01 261 11 Roof02A 91 总数 2212 12 Roof02BGvl 39

表 2 Urban数据的地物真实样本信息 Table 2 The ground truth of samples in each class for Urban dataset

本节利用公开高光谱影像数据集PaviaU和Urban来设计3组实验以验证提出的WPAA用以波段选择的效果。实验采用支持向量机分类器SVM(Support Vector Machine)和K—近邻分类器KNN(K-Nearest Neighbor)来分类波段子集，采用总体分类精度OCA(Overall Classification Accuracy)来定量评价分类效果。SVM分类器采用径向基核函数，其方差和惩罚因子通过交叉验证获得，KNN分类器得邻域大小设定为3。从每个数据集的训练和测试样本中随机抽取10次进行实验，以下是10次独立实验的平均结果。

(1) 综合差异性权重矩阵对WPAA的贡献度分析。本实验对比分析综合差异性权重矩阵W对改进WPAA的贡献度大小。对比方法包括未包含综合差异性权重矩阵的概率原型分析方法PAA(Probabilistic Archetypal Analysis)，基于综合差异性指标的波段选择方法(CDM)和本文提出的WPAA。考虑到波段选择目的为得到高光谱信息量、低波段相关性且类间差异大的波段子集，实验采用3种量化指标包括平均信息熵AIE(Average Information Entropy)、平均相关系数ACC(Average Correlation Coefficient)和平均相对熵ARE(Average Relative Entropy)分别来估计波段子集的光谱信息丰富程度，波段相关关系和波段子集的类间差异性。实验中，PaviaU和Urban数据的波段子集大小k分别人为设定为20和40。

表 3(Table 3)

表 3 不同波段选择方法的定量评价和分类精度对比 Table 3 The contrast in quantitative evaluation and classification accuracy from different band selection methods 数据 方法 定量评价指标 分类器 SVM KNN AIE ACC ARE OCA OCA PaviaU (k=20) PAA 12.643 0.569 1.982 0.918 0.851 CDM 10.541 0.621 1.022 0.894 0.825 WPAA 16.860 0.526 3.254 0.956 0.915 Urban (k=40) PAA 8.682 0.669 10.386 0.952 0.925 CDM 7.841 0.724 6.575 0.915 0.906 WPAA 12.524 0.582 14.067 0.986 0.968

表 3 不同波段选择方法的定量评价和分类精度对比 Table 3 The contrast in quantitative evaluation and classification accuracy from different band selection methods

表3展示PAA、CDM和WPAA在PaviaU和Urban数据集上的定量评价和分类结果。可以看出，PAA方法得到的定量评价指标和分类精度OCA都优于CDM而低于WPAA。3种方法中，WPAA方法在各方面表现最优，能够得到光谱信息量更丰富、波段相关性更小且类间差异性更大的波段子集。这说明通过引入综合差异性权重矩阵到PAA模型中，构建的WPAA模型能够明显提升CDM和PAA的波段选择结果，得到更好的高光谱影像分类结果。因此，综合差异性权重矩阵对保证WPAA的波段选择结果具有重要贡献。

(2) 不同波段数量下的分类结果对比。本实验对比分析不同波段子集大小条件下WPAA和主流对比方法在PaviaU和Urban数据集的分类性能。采用的4种主流对比方法包括基于稀疏的波段选择法SpaBS(Sparse based band selection) (Li和Qi，2011)、稀疏矩阵分解的波段选择法SNMF(Sparse Nonnegative Matrix Factorization) (Li和Qian，2011)、改进稀疏子空间聚类方法ISSC(Improved Sparse Subspace Clustering) (Sun 等，2015b)和稀疏自表达的波段选择方法SSR(Sparse Self Representation) (Elhamifar 等，2012)。

实验中，PaviaU数据集中波段数的选择区间为5—50，Urban数据集中波段数的选择区间为4—45，步长都为5。实验选定每一类地物中20%的真实样本作为训练样本大小来进行分类实验，其余的地物样本作为测试样本。其中PaviaU和Urban数据集中，SpaBS方法的迭代次数设置为5次。利用交叉验证方法，PaviaU数据集中，SNMF方法的正则化因子α和β分别设定为4.0和1.5；Urban数据集中，SNMF方法中正则化因子α和β分别设定为3.5和0.05。

图4展示不同波段子集大小条件下的WPAA和4种对比方法在PaviaU和Urban数据集上的总体分类精度OCA曲线。可以看出，随着波段数量从5增加，5种方法的总体分类精度OCA总体上升。其中，SpaBS的OCA曲线表现最差，较不稳定，明显低于WPAA和其他方法。ISSC方法优于SNMF、SSR和SpaBS，但整体低于WPAA的OCA曲线结果。SNMF和SSR的OCA曲线结果较为接近，都明显优于SpaBS。因此，WPAA在不同波段数量条件下得到的OCA曲线表现最优，分类精度最高。

图 4 不同波段选择方法在PaviaU和Urban数据中的总体分类精度曲线 Figure 4 OCA curves of different band selection methods on PaviaU and Urban datasets

进一步，图5和图6展示WPAA和其他4种方法在选定波段子集数量条件下的分类结果。其中，PaviaU数据中波段子集k的大小为20，Urban数据集中波段子集k的大小为40。可以看出，WPAA的分类图结果优于其他4种方法。

图 5 不同波段选择方法的PaviaU数据的SVM分类结果 Figure 5 SVM classification maps of different band selection methods on PaviaU dataset

图 6 不同波段选择方法的Urban数据的SVM分类结果 Figure 6 SVM classification maps of different band selection methods on Urban dataset

(3) 不同训练样本大小下的分类结果对比。本实验对比分析不同训练样本条件下WPAA和其他4种对比方法(SpaBS、SNMF、ISSC和SSR)的分类结果。其中，PaviaU和Urban数据的训练样本随机采样于地物真实样本集，采样比率的选择区间为[0.05, 0.08, 0.1, 0.13, 0.15, 0.18, 0.2, 0.25, 0.3, 0.4, 0.5]。PaviaU数据中，波段数量k设定为20；Urban数据中，波段数量k设定为40。IKSSR和其他4种波段选择方法的其他参数设置都与实验(1)保持一致。

图 7 PaviaU和rban数据中不同训练样本大小下的各种波段选择方法的总体分类精度OCA曲线 Figure 7 OCA curves of different band selection methods with different size of training samples per class on PaviaU and Urban datasets

图7为不同地物训练样本大小条件下的5种方法在PaviaU和Urban数据中得到的总体分类精度OCA曲线结果。考虑到SVM和KNN分类结果趋势的一致性，图7仅列出SVM的分类结果。可以看出，随着单一地物训练样本大小的增加，各种波段选择方法得到的分类结果提升较为明显。对于任一高光谱数据集，WPAA得到的总体分类精度OCA始终高于其他4种方法包括SpaBS、SNMF、ISSC和SSR。类似于实验(1)，SpaBS得到的总体分类精度OCA明显低于其他方法。SNMF和SSR的分类精度较为接近，但两者的OCA曲线低于ISSC的OCA结果。ISSC的OCA曲线高于SpaBS、SNMF和SSR，但低于WPAA方法。因此，WPAA在不同训练样本条件下得到的总体分类精度OCA最高，优于其他4种对比波段选择方法。

本文提出加权概率原型分析WPAA方法来解决高光谱影像的波段选择问题。WPAA方法能够考虑波段之间的差异性和高光谱遥感成像的量子特性，引入贝叶斯框架理论和波段权重矩阵来改进传统的原型分析模型。整合WPAA方法采用迭代优化的策略来转换非凸优化问题，运用交替方向乘积方法来求解得到局部最优的稀疏系数矩阵来得到合适的波段子集。利用PaviaU和Urban两个公开高光谱数据集，对比SpaBS、SNMF、ISSC和SSR 4种波段选择方法来验证WPAA的波段选择效果。实验结果表明，引入的综合差异性权重矩阵对保证WPAA的波段选择结果具有非常重要的作用。另一方面，不同波段子集和不同训练样本下的WPAA的OCA曲线都高于4种主流方法，展现其波段子集用于高光谱影像的更高精度分类的能力。因此，WPAA能够选择适合高光谱影像的波段子集，更好满足分类的应用需求。然而，本文的研究结果需要下一步的工作来进行完善和优化。在后续的工作中将考虑采用泊松分布来模拟噪声的概率分布来优化WPAA的波段选择结果。同时，将考虑采用波段高维空间分布的图谱理论来改进WPAA模型，提升WPAA的波段选择效果来满足实际高光谱影像数据处理需求。