基于StOMP稀疏方法的高光谱图像目标检测

赵春晖，靖晓昊，李威

(哈尔滨工程大学信息与通信工程学院，黑龙江哈尔滨150001)

高光谱遥感技术是一门综合性的高新技术，在理论、技术和应用上都得到了长足的发展，高光谱遥感同时利用空间图像和光谱特征来获取地物信息，实现了图谱合一［1］。因为高光谱图像同时包含空间信息与光谱信息，所以其在军事目标检测、运动目标识别等方面有着特有的优势。目标检测可以被视为将图像中的像元标记为目标或是背景。随着遥感技术的发展，高光谱图像目标检测技术在农业、军事等领域发挥着越来越重要的作用［2-3］。目前比较成熟的检测算法［4-5］包括异常目标检测(RX 算法)［6］、支持向量机(support vector machines，SVM)［7］、光谱匹配滤波(spectral matched filter，SMF)［8］、匹配子空间检测(matched subspace detector，MSD)［9］和自适应子空间探测(adaptive subspace detector，ASD)［10］等。近年来，稀疏表示方法得到了人们越来越广泛的关注，也逐渐开始有学者将稀疏表示的思想引入到高光谱图像目标检测领域，为高光谱图像目标检测开辟了一条新的思路。随着稀疏方法的应用，新的问题也渐渐浮现出来:稀疏表示方法中大量的矩阵运算有着很高的计算成本，而高光谱图像的数据量与一般的二维全彩图像相比又要大出很多，这就使得高光谱图像的稀疏目标检测方法的运算量非常巨大，数据处理时间往往很长，也就导致了稀疏方法在某些高光谱图像目标检测应用领域的局限性。基于以上原因，本文提出了基于StOMP重构算法［11］的高光谱图像目标稀疏检测方法，相比于传统算法，StOMP算法减少了大量迭代次数，大大降低了运算成本，提升了检测速度;同时，在保证检测精度不明显下降前提下，在某些复杂场合能够进一步提高检测效果。

1 高光谱图像稀疏表示与目标检测

稀疏表示的基本思想就是将原始信号X表示为字典D与系数A乘积的形式:其中，要求字典D是过完备的，即K≫N，K为字典D中原子的个数，N为信号的长度。由于算法并不要求字典中原子的正交性，所以系数A并不唯一，一般要找出最为稀疏的一组系数，这样便可以通过少数几个大系数来展现信号的本质特征。在实际应用中，可以用几个系数来很好的描述一个信号，简化后续处理工作量。

在高光谱图像中，任何一个像元的光谱信息均可以表示为大量端元中的某些原子光谱的线性组合;此外，高光谱图像由于包含了丰富的光谱信息，数据量相较于一般的全彩图像要大出很多。不难看出，高光谱数据非常适合进行稀疏分解处理，而事实上，高光谱图像分解后得到的稀疏系数也往往具有很高的稀疏度。

令X是一组含有B个波段的高光谱图像数据，D= di{ }i=1，…，N是含有N个原子的过完备字典。由上面的叙述可知，高光谱图像中的某个B维像元x可以表示为字典D中原子的线性组合:

式中:过完备字典D表示B×N矩阵;a为一个未知的N维向量，表示在上述线性组合中每个原子的权重。根据稀疏表示理论，a是稀疏的，即a中只包含少数非零元素。本文中字典包含的原子来自于高光谱图像中随机选取的像元。

高光谱图像目标检测可以被视为一种特殊的二元分类，即图像中只有背景类和目标类。因为图像中的像元光谱可以由字典原子的线性组合表示，可以根据表示像元光谱所使用的字典原子属于背景子字典还是目标子字典来判断该像元的归属［12］。对于任意的待检测像元 x，可以假设它坐落于一个由背景子空间和目标子空间共同张成的空间里。它的光谱则可以近似地表示为背景子字典Db和目标子字典Dt中对应的训练样本的线性组合:

式中:字典D是由目标子字典Db和背景子字典Dt构成的B× Nb+Nt

( )矩阵;a是由ab和at组合而成的Nb+Nt维权重系数。

在稀疏模型中，不需要假设目标和背景的分布;又由于目标像元和背景像元的光谱特性存在区别，使其分布于不同的子空间中。稀疏向量a由背景权重系数ab和目标权重系数at组合而成，若x是一个目标像元，则ab是零向量，at是一个稀疏向量;若x是一个背景像元，则ab是一个稀疏向量，at是零向量。因此，根据像元x稀疏表示系数a的非零元位置就可以判别该像元是背景像元还是目标像元［13］。

求解像元光谱x的稀疏表示的系数向量a就是求解下式的最优化问题:

式中:‖·‖0代表求解l0范数，即非零元素的个数，但是l0范数最小化问题是一个NP-hard问题。由于向量a稀疏的特性，可以用求解l1范数的最小值问题来逼近最小化l0范数的问题，近似后如下

由上式求得的a是在字典D下分解像元光谱x得到的最稀疏权重系数，它由背景权重系数向量ab和目标权重系数向量at组合而成。根据式(5)由重构算法求解出a后，便可以通过比较重构残差:

的值来判断x的归属类别。因此，将检测器的输出设为

给定门限δ，如果R x()＞δ，则判定x为目标像元，否则，则判定像元x为背景像元。

2 StOMP目标检测

准确地求解稀疏系数a是高光谱图像稀疏目标检测算法的关键环节之一，常见的重构算法有匹配追踪(matching pursuits，MP)算法、正交匹配追踪算法(orthogonal matching pursuit，OMP)等。但由于MP算法的字典中原子数目太大，而且每次匹配只选取一个最佳原子，就使得整个选取过程中匹配次数过多，继而增加算法的计算复杂度;而OMP算法方法利用Gram-Schmidt正交化过程将投影方向正交化，从而提高了追踪的效率［14］，但是其正交化过程引入了新的计算开销，特别是对于高光谱图像信号，计算量仍然巨大。

StOMP算法将OMP算法进行进一步的简化，大大提高了运算速度。这一算法的本质在于每次匹配追踪时选出的是多个匹配原子而不是单个原子，减少了匹配次数［7］。算法的流程如下:

第1步:令m=0，下标集合Im=φ，首先将信号f投影到D一个向量gr∈D上，计算出余项Rf:

因Rf与gr正交，故有

选取匹配原子，引入阈值极限因子t0，并找到满足条件原子的下标集合:

利用Gram-Schmidt正交化过程将这组原子正交化，并定义正交化后的向量为

将f投影到uii∈I0( )上，得到

在找到这些匹配原子后，与上一次迭代得到了原子下标集合合并，从而更新所选的原子下标集合Im=Im∪Im-1在 ui:i∈Im-1{

}的基础上对新引入的原子继续正交化。之后仿照第一步的过程得到

迭代结束后，设共迭代了M次，则此时有

StOMP稀疏目标检测算法流程如图1所示。

图1StOMP算法流程图Fig.1 StOMP algorithm flowchart

3 仿真结果及分析

3.1 高光谱图像实验数据

为了验证本文所提算法的有效性，共使用2幅高光谱图像数据进行仿真实验，这2幅图像数据的第50波段示意图，如图2所示。

图2 仿真实验数据(第50波段)Fig.2 Simulation experiment dataset(band 50)

图像数据图2(a)是来自文献［4］的一幅合成数据。该图像的尺寸为30×30，包含6个目标，目标尺寸分别为 3×3、4×4 和 5×5，目标排列如图 2(a)所示。这幅图像包含了126个波段的高光谱数据。这幅合成图像的构造过程如下:先在本文所使用的第3幅飞机场图像数据的所有目标像元中随机的选择了100个像元，按照图2(a)中的目标位置进行排列，再在飞机场图像数据的所有背景像元中随机的选择800个像元排列在目标像元周围，构成了这一30×30的共900个像元的合成图像数据。

图像数据图2(b)是利用先进的机载可见光/红外成像光谱仪(airborne visible infrared imaging spectrometer，AVIRIS)采集的数据。AVIRIS是采用推扫成像方式的成像光谱仪，在 0.4～2.45 μm 的波长范围内获取224个波长处的空间图像信息，波长间隔为10 nm［15］。本实验所用图像是美国圣地亚哥机场的一部分，它覆盖了从可见光到近红外的光谱范围，去除水的吸收带和信噪比较低的波段后，余下的126个波段参与仿真实验［4］。所用实验图像大小为100×100，图中包含了38个待检测目标，如图2(b)所示。

本文仿真环境:处理器:64位双核Intel i5-2450 M，主频 2.50 GHz;内存 4 G;操作系统:Microsoft Windows7;程序运行平台:MATLAB 2011b。

3.2 实验结果与分析

为了验证本文算法的效果，选择了MP算法、OMP算法与本文所采用的StOMP算法作比较，分别将这3种算法应用于高光谱图像目标稀疏检测模型。

首先，将3种检测模型应用于合成图像2(a)，3种算法的检测效果如图3所示。

图3 合成数据仿真检测结果Fig.3 Detection results of manmade dataset

从图3可以看出，对于合成数据，由于数据相对简单，OMP和StOMP2种算法都可以准确地检测出目标像元的位置，而MP算法的检测效果明显劣于另外2种算法。

表1给出了3种算法在处理合成图像时所消耗的时间。从表中可以看出，与MP和OMP算法相比，由于StOMP算法在匹配的过程中每次取出多个而不是一个原子，减少了迭代次数，故StOMP算法消耗的时间比MP和OMP算法都要少，运算速度有明显优势。

表1 2种图像运算时间比较Table 1 Computing time of two images

为了进一步论证算法的效果，再将检测模型应用于真实飞机场图像2(b)，3种算法的检测效果如图4所示。

从图4可以看出，对比合成图像，3种算法应用于真实图像的检测精度差异更为明显。虽然StOMP算法的迭代次数相较于OMP算法减少了很多，但由于匹配时每次取出多个原子，追踪时匹配过的原子总数并不会减少很多，所以OMP算法与StOMP算法的检测效果相仿，均能正确地检测出图像中的每个目标，而MP算法由于无法保证追踪过程中的正交性，检测效果最差。

图4 飞机场数据仿真检测结果Fig.4 Detection results of real airport dataset

表1给出了3种算法在处理真实飞机场图像时所消耗的时间。从表1可以清楚地的看到，StOMP算法消耗的时间远远少于MP与OMP算法，对比表1发现，越是大型的数据，StOMP算法速度快的优势就越是明显。这是因为，真实图像的尺寸要比合成图像大很多，矩阵的尺寸越大，运算量就越大;而匹配过程中，每次迭代都要进行大量的矩阵乘法、矩阵求逆等运算，所以每减少一次迭代，对减少运算量做出的贡献都是很大的。

利用接收机工作特性(ROC)曲线客观的评价3种算法的检测性能。接收机工作特性曲线用于描述不同检测阈值下检测概率pd与虚警概率pf之间的变化关系，提供对算法检测性能的定量分析。检测概率pd为检测到的真实目标像素数目与地面真实目标像素数目的比值;虚警概率pf为检测到的虚警像素数目同整幅图像像素数目总和的比值。3种算法针对真实飞机场图像数据的ROC曲线如图5所示。

图5 飞机场数据接收机工作特性曲线Fig.5 ROC curves of airport dataset

从ROC曲线可以看出，StOMP算法与OMP算法的检测性能相近，而在阈值小于0.1时，StOMP算法的效果要略优于OMP算法，这就意味着StOMP算法对于虚警概率要求较高的情况有着更好的适用性，此外二者的检测效果均大大优于MP算法。

4 结束语

本文介绍了高光谱图像目标稀疏检测理论，并分析了其存在的不足，针对计算量过大的问题，提出了利用StOMP算法进行目标检测的方法，并取得了较好的检测效果。应用StOMP算法可以在保证检测精度不下降的前提下，大幅提升检测速度，在检测速度要求比较高的场合具有广阔的研究应用价值。

［1］王立国，赵春晖.高光谱图像处理技术［M］.北京:国防工业出版社，2013.

［2］MATTEOLI S，DIANI M，CORSINI G.A tutorial overview of anomaly detection in hyperspectral images［J］.IEEE A＆E Systems Magazine，2010，25(7):5-27.

［3］DATT B，MCVICAR T R，Van NIEL T G，et al.Preprocessing EO-1 Hyperion hyperspectral data to support the application of agricultural indexes［J］.IEEE T Geosci Remote，2003，41(6):1246-1259.

［4］ZHAO Chunhui，LI Xiaohui，REN Jinchang，et al.Improved sparse representation using adaptive spatial support for effective target detection in hyperspectral imagery［J］.International Journal of Remote Sensing，2013，34(24):8669-8684.

［5］BANERJEE A，BURLINA P，DIEHL C.A support vector method for anomaly detection in hyperspectral imagery［J］.IEEE T GEOSCI REMOTE，2006，44(8):2282-2291.

［6］KRAUT S，SCHARF L L.Adaptive subspace detectors［J］.IEEE T Signal Proces，2001，49(1):1-16.

［7］赵春晖，成宝芝，杨伟超.利用约束非负矩阵分解的高光谱解混算法［J］.哈尔滨工程大学学报，2012，33(3):377-382.

ZHAO Chunhui，CHENG Baozhi，YANG Weichao.Algorithm for hyperspectral unmixing using constrained nonnegative matrix factorization［J］.Journal of Harbin Engineering University，2012，3(33):377-382.

［8］MANOLAKIS D，SHAW G.Detection algorithms for hyperspectral imaging applications［J］.IEEE Signal Processing Magazine，2002，19(1):29-43.

［9］赵春晖，李杰，梅峰.核加权RX高光谱图像异常检测算法［J］.红外与毫米波学报，2010，29(5):378-382.

ZHAO Chunhui，LI Jie，MEI Feng.A kernel weighted RX algorithm for anomaly detection in hyperspectral imagery［J］.J Infrared Millim Waves，2010，29(5):378-382.

［10］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-3985.

［11］DONOHO D L，TSAIG Y，DRORI I，et al.Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit［J］.IEEE Transactions on Information Theory，2012，58(2):1094-1121.

［12］宋相法，焦李成.基于稀疏表示及光谱信息的高光谱遥感图像分类［J］.电子与信息学报，2012，34(2):268-272.

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-272.

［13］CHEN Y，NASRABADI N M，TRAN T D.Sparse representation for target detection in hyperspectral imagery［J］.IEEE J Sel Top Quant，2011，5(3):629-634.

［14］PAT Y C，REZAIIFA R，KRISHNAPRAS P S.Orthogonal matching pursuit:recursive function approximation with applications to wavelet decomposition［C］//Proceedings of the 27th Annual Asilomar Conference on Signals，Systems，and Computers， Signals， Systems and Computers.Los Alamitos，USA，1993:40-44.

［15］赵春晖，李晓慧，朱海峰.空间4-邻域稀疏表示的高光谱图像目标检测［J］.哈尔滨工程大学学报，2013，34(9):1-8.

ZHAO Chunhui，LI Xiaohui，ZHU Haifeng.Hyperspectral imaging target detection algorithm based on spatial 4 neighborhood for sparse representation ［J］.Journal of Harbin Engineering University，2013，34(9):1-8.