基于随机矩阵的高光谱影像非负稀疏表达分类

孙伟伟，刘 春，2，施蓓琦，李巍岳

（1.同济大学 测绘与地理信息学院，上海200092；2.现代工程测量国家测绘地理信息局重点实验室，上海200092）

高光谱影像通过成像光谱仪采集地物的光谱特征，具有波段众多且光谱分辨率高的特点.通过对高光谱像元光谱特征分类，能够识别地物间的细微差异.高光谱影像分类目前广泛应用于环境监测［1］、土地利用分析［2］及资源勘查［3］等方面.

近几年来，压缩感知（compressive sensing，CS）的出现和兴起，为高光谱影像的分类研究提供了新的研究思路.CS理论指出，如果一个信号在某个变化上稀疏，那么可以用一个与变换基不相关的测量矩阵将稀疏的高维信号投影到低维空间上，通过求解一个优化问题高概率地还原原始信号［4］.在高光谱研究领域，当前学者已经在CS用于分类方面做出一些探索.Haq等［5］提出基于稀疏表达和L1范数最小的高光谱影像分类方法.随后，Chen 等［6］提出基于光谱稀疏表达和重构的高光谱影像分类模型，并引入窗口邻域及K－邻域来描述空间信息以提高分类精度.Castrodad等［7］研究附加像元的空间邻域信息至高光谱影像稀疏表达分类，并扩展稀疏表达分类法至非监督领域.近期，国内学者利用字典学习构造像元的稀疏表达并结合随机森林法来实现影像分类［8］；通过光谱维和空间维的联合表达和约束构建每类训练样本的稀疏字典及稀疏表达，并结合最小重构误差和邻域相关性约束求解系数向量来实现分类［9］.然而，高光谱影像稀疏表达分类模型仍存在两个主要问题：①测量矩阵是所有类别训练样本或纯净像元的光谱向量排列矩阵，列向量间尤其同类地物间相关性强，不满足CS 要求的限制等距特性（restricted isometry property，RIP）准则；②分类模型未考虑光谱向量在高维样本空间中重构系数向量的非负性.以上两个问题关系高光谱影像稀疏表达分类模型理论的严密性和准确性，严重影响高光谱影像的分类精度.

本文提出随机矩阵－非负稀疏表达（random matrix－nonnegative sparse representation， RMNSR）分类模型，通过引入随机映射来构造随机测量矩阵以提高原始测量矩阵更好满足RIP 条件，同时结合限定重构系数向量的非负性，最终提高CS理论用于高光谱影像分类的精度.同时，通过两个高光谱数据实验，在三种常用的重构算法下，对比改进RMNSR 模 型 与 常 规 的 稀 疏 表 达 （sparse representation，SR）分类模型对分类的精度差异，证明所提模型的有效性.

1 压缩感知理论

CS理论认为，如果一个信号具有稀疏或可压缩的特性，那么可以采用线性和非自适应的测量过程，从较少的测量值中准确地恢复出原始信号［10］.假设将一组M维不含噪声的信号x∈RM×1在某组基（φi为M维列向量）下展开，表示如下：

其中θi＝〈x，φi〉＝φTix，Φ＝［φ1，φ2，…，φM］∈RM×M为稀疏矩阵且满足正交条件ΦΦT＝ΦTΦ＝I.系数向量θ＝［θ1，θ2，…，θM］T为K－稀疏的，即非零系数的个数K≪M.在此基础上，采用一个与Φ 不相关的测量矩阵ΨN×M（N≪M），通过式（2）对x 实现压缩观测：

可以得到N个线性观测y∈RN×1，并包含重构信号x的信息，其中N为观测值的个数.将式（1）代入式（2），可以得到：

其中A＝ΨΦ 为信息算子.式（3）是一个病态方程，但只要矩阵A 中任意2K列线性独立，那么至少存在一个K－稀疏的系数向量θ 满足y＝Aθ［11］.因此，通过求解一个非线性优化问题（4）：

可以获得稀疏系数向量θ 的唯一解，进而高精度地重建信号x.式（4）的非凸优化问题通常采用重构算法 来 求 解，如L0 范 数 法［12］、基 追 踪 法［13］（basis pursuit，BP）和 正 交 匹 配 追 踪［14］（orthogonal matching pursuit，OMP）.重构方法具有各自不同的特点，其选择依赖于非线性优化目标方程.

2 基于RM－NSR的高光谱影像分类

假设归一化的高光谱影像数据集为X＝｛xi｝∈RD×N，其中D为像元的个数，N为波段数.假设X共包含c个类别，训练样本总数为M，且第r个类别包含的训练样本为，其中mr为第r类别的训练样本个数.高光谱影像的RM－NSR 分类目的是基于随机测量矩阵，通过CS理论来重构稀疏的非负系数向量，最终确定测试样本中每个像元x的类别标签.

2.1 光谱向量在样本空间的非负稀疏表达

光谱向量稀疏表达分类模型认为高光谱影像中同类地物的像元都位于统一的特征空间.因此如果像元x位于第r类，则其光谱向量可看作位于由第r类的训练样本张成的子空间，即x 可表达为训练样本构成的矩阵的线性组合，

其中Ψr为第r类训练样本构成的矩阵，每一列为该类中训练样本的光谱向量，θr中元素为第r类训练样本中光谱向量的权重系数.如果将x 在所有类别的训练样本构成的高维空间中展开，则

式中：Ψ 为所有类别的训练样本构成的测量矩阵；Θ为像元x在训练样本中所有光谱向量所对应的权重系数向量；M为权重系数的个数，等于训练样本的总数；N为光谱向量的维数.在理想的情况下，如果位于第r类的像元x 可以由其位于同一类别的光谱向量线性表达，则Θ 为稀疏向量考虑光谱响应值及权重系数的非负性，限定第r类训练样本对应的权重系数＝1，2，…，mr，这不同于常规SR 模型.

2.2 随机测量矩阵

CS理论理论指出，当测量矩阵Ψ 满足RIP 特性时，才能用较少的测量值x 严格恢复原始的高维稀疏系数向量Θ.通常，可以用Ψ 的不自相关性来检验其是否满足RIP特性.而式（6）中，Ψ 为所有训练样本的光谱向量排列矩阵，每一列光谱向量及每一维的光谱值间相关性很强，因此并不满足RIP准则.理论证明［10］，如果测量矩阵中每个元素服从均值为0且方差为1／N的严格子高斯分布，该测量矩阵满足RIP准则；当时可以高概率地重构原始信号，其中N为测量值个数，M为高维信号的维数，K为稀疏度，κ1为任意值.因此，采用随机矩阵来逼近满足严格子高斯分布的测量矩阵.引入随机矩阵RP×N到式（6）两侧，得到，

其中x′为随机映射后的低维光谱向量，维数为P；Ψ′为随机测量矩阵，P≪N，N为原始光谱向量的维数.矩阵R 满足严格子高斯分布，应用随机矩阵到式（6）的两侧可降低Ψ 中每一列的相关性，并使得测量矩阵Ψ′的每一列满足严格子高斯分布，进而改善测量矩阵满足RIP条件.同时，随机矩阵通过投影能够减少光谱向量x的维数，降低重构系数向量的计算量.常用的随机映射矩阵有均值为0且方差为1的高斯分布矩阵或伯努利随机矩阵.此外，可以看出投影维数P对式（7）中Θ 及分类结果产生影响.

2.3 系数向量的非负重构及高光谱影像分类

式（7）中，由于P≪N≪M，系数向量Θ 并不存在唯一解，其求解过程是一个非凸问题.通常，Θ 的非零元素通过求解式（8）的非线性优化问题得到：

式（8）可通过平滑L0范数的方法来求解，利用稳定的连续方程来逼近不连续的L0范数，通过最小化连续方程得到最优L0范数［12］.

另一方面，通常采用L1范数代替L0范数来松弛式（8）得到凸优化问题，同时考虑到实际优化过程中的重构误差，得到：

其中，σ为设定的重构误差，通过交叉验证获得.对式（9）的求解称为基追踪（L1范数最小）［13］，通常采用线性凸规划保证L1范数收敛到全局最小，得到最优解.此外，限定系数向量的稀疏度K，式（8）的L0范数最小也可转换为

式（10）可通过贪婪算法，如正交匹配追踪［14］，通过迭代选择局部最优解来逼近原始信号，得到重构理想情况下，根据恢复的系数向量中的非零元素位置，可判别该像元所属的地物类别.然而由于噪声的影响，非零元素可能出现在多种地物类别对应的位置上.因此，采用每一类训练样本对应的重构误差来判断该地物类别［15］.通过式（11）中二阶范数计算得到该像元对每一类训练样本的重构误差，其最小值对应的类别则为该像元的地物类别，

其中：r＝1，2，…，c.

2.4 高光谱影像RM－NSR的分类流程

高光谱影像RM－NSR分类的流程如下：①基于归一化的高光谱训练样本，利用式（6）构造矩阵Ψ；②构造随机矩阵R，利用式（7）得到随机测量矩阵Ψ′及随机矩阵投影后的低维光谱向量x′；③针对式（8）—（10）中不同的目标优化问题，采用对应的算法得到恢复的高维稀疏系数向量；④对每一像元，计算得到每一类别的重构误差，利用式（11）中的最小重构误差确定其类别.

3 实验和分析

3.1 实验数据

Urban高光谱数据是从美国陆军地理空间中心获取的HYDICE 数据.数据采集于1995年10 月，空间分辨率为2m，光谱分辨率为10nm.影像大小为307×307像素，覆盖美国德克萨斯州科帕拉斯区域（靠近胡德堡），如图1所示.对原始的210波段数据进行预处理，移除噪声波段区间［1－4，76，87，101－111，136－153，198－210］，剩余162个波段，包含22种主要地物.PaviaU 高光谱数据来自于西班牙巴斯克大学计算智能课题组，影像覆盖帕维亚大学区域，共103波段，空间分辨率为1.3m，如图2所示.影像为较大数据集中的一部分，包含350×340像素，波段数为103，包含9类地物（包括阴影）.数据1和数据2中每一类地物中训练和测试样本的信息见表1.

表1 Urban和PaviaU 数据中训练和测试样本信息Tab.1 The information of training and testing samples in Urban and PaviaU datasets

3.2 实验分析

采用RM－NSR 模型对Urban数据和PaviaU 数据进行分类实验，验证RM－NSR 模型提高分类精度的有效性.同时，研究随机投影的维数P对分类精度的影响，以更好把握分类模型.为了综合评价RMNSR 模型的分类结果，求解模型下的三种常用的目标优化问题，得到非负系数向量.其中，式（7）中采用均值为0方差为1／P的随机矩阵，式（8）采用非负约束的平滑L0重构算法（R＋L0），式（9）采用非负约束的基追踪算法（R＋L1），式（10）采用非负约束的正交投影追踪算法（R＋OMP）.在此基础上，对比传统SR 模型采用的L0、基追踪（L1）和OMP 算法得到的分类结果.

（1）分类结果分析

Urban数据的RM－NSR 模型和传统SR 模型的分类结果对比见表2.其中，RM－NSR 模型中投影维数P为60.R＋L0 和L0 中最小重构误差阈值为0.000 1，迭代次数为100，重构误差下降因子为0.5；R＋L1和L1中重构误差阈值为0.000 1；R＋OMP和OMP中，重构误差阈值为0.000 1，迭代次数为10.RM－NSR 的R＋L0比SR 模型中L0 的平均分类精度高5.29%，R＋L1比L1分类精度高5.06%，而R＋OMP比OMP的分类精度高5.08%.单一地物分类结果中，大多数地物的R＋L0分类精度都比L0都较大幅度提高，如第4类的Vegpasture和第6类的VegTrees01.每一类地物的R＋L1的分类结果比L1也有很大改进，尤其第3类的Concrete01.此外，R＋OMP的单个地物分类精度很多达到100%，相比OMP有很大提升.三种不同重构方法的分类结果对比说明，RM－NSR 模型能够明显改善传统SR模型对Urban 数据的分类精度，平均提高约5.14%.

表2 Urban数据分类结果对比Tab.2 The contrast in classification results for Urban dataset

PaviaU 数 据 的RM－NSR 模 型 与 传 统SR 模 型的分类结果对比见表3.其中，RM－NSR 模型中随机矩阵的投影维数P为50.R＋L0和L0中最小重构误差阈值为0.001，迭代次数为80，重构误差下降因子为2；R＋L1和L1中重构误差阈值为0.001；R＋OMP和OMP中，重构误差与之为0.001，迭代次数为15.由表3 可以看出，RM－NSR 的R＋L0 比L0的平均分类精度提高3.87%，R＋OMP 比OMP 的分类结果精度高4.33%，而R＋L1相比L1分类精度增幅最大，为5.64%.第1 类Asphalt和第3 类Gravel的R＋L0分类精度相比L0增幅最大，分别为11.39%和10.73%.R＋L1和L1的每一类的分类结果对比中，第5类Painted metal sheets的分类结果提升明显，大于22%.相比R＋L0和R＋L1，R＋OMP的单类精度增幅较为均匀.三类重构方法的分类结果对比说明，相比SR 模型，RM－NSR 能够明显提高PaviaU 数据的分类精度，平均约4.61%.

表3 PaviaU 数据分类结果对比Tab.3 The contrast in classification results for PaviaU dataset

（2）随机矩阵中投影维数P的影响

实验证明随机矩阵R能够明显改善传统SR 模型的分类结果，然而式（7）中随机矩阵的投影维数P对结果有很大影响.因此，需要研究P对分类精度的关系，确定合适的P来保证分类结果的改善.图3为Urban和PaviaU 数据的RM－NSR 模型中投影维数P与平均精度的关系图.其中，R＋L0、R＋L1和R＋OMP的重构误差及迭代次数设置与前面实验保持一致，P的区间为［3，100］.

图3 中随机矩阵的投影维数P 与平均分类精度的关系Fig.3 The relationship between the projected dimension P and the average classification accuracy for Urban and PaviaU datasets

从图3a看出，当投影维数P较小（小于3）时，三种重构方法得到的分类精度都很低.随着P增大，平均分类精度开始急速提升，在P＝22左右达到较高值且大于SR 模型的分类精度.然而，在P＞22后，分类精度增长开始缓慢，随着投影维数增加而波动性地增长，呈现较为稳定的趋势.同样，图3b 的PaviaU 数据中，当P较小时，分类精度随着投影维数的增加而急剧增大，在P＝10（类别个数为9）左右达到较高值；在P＞10后，随着投影维数增大，分类精度较为稳定并伴随小幅波动，总体大于SR 模型的分类结果.总结得到以下结论：过小的投影维数将产生较低的分类精度；当投影维数P明显小于类别数时，随着P的增大，分类精度急剧增大并在一定数值后达到较高的值；此后，随着P的增大，分类精度整体缓慢上升同时呈现小幅的波动趋势.考虑实际RM－NSR 模型对高光谱影像分类精度提高的要求，应选取较大的且大于类别个数的投影维数.

4 结论和展望

本文提出随机矩阵－非负稀疏表达（RM－NSR）分类模型来改进常规的稀疏表达（SR）分类模型，以改善高光谱影像的分类精度.通过采用随机矩阵来改善各类别的光谱向量的排列矩阵，形成新的每列满足严格子高斯分布的测量矩阵，从而改善测量矩阵满足的RIP条件.同时，对重构的高维稀疏向量附加非负约束，提高系数向量的可解释性以符合实际.基于两个不同的高光谱数据，针对三种常用的优化问题，对比RM－NSR 和SR 模型的三种重构方法得到的分类结果.实验表明，对三种重构方法，RMNSR 模型能够明显提升SR 模型用于高光谱影像分类的平均精度.此外，随机矩阵的投影维数对平均精度的影响研究表明，较小的投影维数将导致较低的分类精度，平均精度随着投影维数的增加先急剧增大后增长缓慢.在实际应用中，应选用较大的投影维数来保证RM－NSR 模型用以提高分类精度.然而，RM－NSR 分类模型中的重构误差的设定通过人为设定来获得.在接下来工作中，将研究重构误差σ的优化设置来进一步推广本文的方法.

致谢：感谢美国马里兰大学帕克分校数学系诺伯特维纳中心的J J Benedetto教授、Wojciech Czaja副教授、王蓉蓉和Yen－ming Mark Lai博士对本论文的指导和帮助.

［1］ Lee B S，McGwire K C，Fritsen C H.Identification and quantification of aquatic vegetation with hyperspectral remote sensing in western Nevada rivers，USA［J］.International Journal of Remote Sensing，2011，32（24）：9093.

［2］ Yang H，Ma B，Du Q，et al.Improving urban land use and land cover classification from high－spatial－resolution hyperspectral imagery using contextual information［J］.Journal of Applied Remote Sensing，2010，4（1）：041890.

［3］ Bishop C A，Liu J G，Mason P J.Hyperspectral remote sensing for mineral exploration in Pulang，Yunnan Province，China［J］.International Journal of Remote Sensing，2011，32（9）：2409.

［4］ Donoho D L.Compressed sensing［J］.IEEE Transactions on Information Theory，2006，52（4）：1289.

［5］ Haq Q，Shi L，Tao L，et al.A L1－minimization based approach for hyperspectral data classification［C］／／Proceedings of 2010 International Conference on Remote Sensing.Hangzhou：［s.n.］，2010：139－142.

［6］ Chen Y，Nasrabadi N M，Tran T D.Hyperspectral image classification using dictionary－based sparse representation［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49（10）：3973.

［7］ Castrodad A，Xing Z，Greer J B，et al.Learning discriminative sparse representations for modeling，source separation，and mapping of hyperspectral imagery［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49（11）：4263.

［8］ 宋相法，焦李成.基于稀疏表示及光谱信息的高光谱遥感图像分类［J］.电子与信息学报，2012，34（2）：268.SONG Xiangfa，JIAO Licheng.Classification of hyperspectral remote sensing image based on sparse representation and spectral information［J］.Journal of Electronics ＆Information Technology，2012，34（2）：268.

［9］ 宋琳，程咏梅，赵永强.基于稀疏表示模型和自回归模型的高光谱分类［J］.光学学报，2012，32（3）：314.SONG lin，CHENG yongmei，ZHAO yongqiang.Hyperspectrum classification based on sparse representation model and autoregressive model［J］.ActaOpticaSinica，2012，32（3）：314.

［10］ Davenport M A，Duarte M F.Introduction to compressed sensing［J］.Electrical Engineering，2011，93：1.

［11］ Donoho D L，Elad M.Optimally sparse representation in general（non－orthogonal）dictionaries via l1 minimization［R］.San Francisco：Department of Statistics of Stanford University，2002.

［12］ Mohimani H，Babaie－Zadeh M，Jutten C.A fast approach for overcomplete sparse decomposition based on smoothed norm［J］.IEEE Transactions on Signal Processing，2009，57（1）：289.

［13］ Yang J，Zhang Y.Alternating direction algorithms for L1－problems in compressive sensing［J］.SIAM Journal on Scientific Computing，2011，33（1）：250.

［14］ Tropp J A，Gilbert A C.Signal recovery from random measurements via orthogonal matching pursuit［J］.IEEE Transactions on Information Theory，2007，53（12）：4655.

［15］ Yang A Y，Wright J，Ma Y，et al.Feature selection in face recognition： a sparse representation perspective ［R］.Berkeley：University of California，Berkeley，2007.