粒子群优化的高光谱图像快速稀疏分解算法∗

王 丽 孙长杰王 威

（1.西安航空学院 西安 710077）（2.中国飞行试验研究院 西安 710089）

1 引言

高光谱图像（Hyperspectral Images，HSIs）包含丰富的空间几何信息和光谱特征信息，适用于战场地形侦察［1］、目标探测与识别［2～3］、分类［4～5］等军事和民用领域。高光谱成像光谱仪的研制［6］、高光谱图像解译研究［7］，已成为备受关注的前沿热点。应用领域的不断扩展，要求高光谱图像能够提供更细致的目标信息。问题是，空间和光谱分辨率越高，对系统的数据存储和传输能力的挑战更高。将信号稀疏分解理论用于高光谱图像［8］，便于对高光谱图像进行后续处理。高光谱图像波段多，如何寻找计算效率高的稀疏分解算法成为研究的关键。

经典的稀疏分解算法是正交匹配追踪（Orthog-onal Matching Pursuit，OMP）［9］算法，该算法能够从原子库中依次查找到能够对信号进行稀疏表示的原子，重构精度高。高光谱图像特征复杂，采用冗余度高的冗余字典才能完成其稀疏表示。正交匹配追踪算法在处理此类问题时，计算时间在现有计算条件下无法忍受。鉴于进化算法［10］寻优能力强的特点，本文采用粒子群算法（Particle Swarm Optimization，PSO）［11］和Hermitian求逆引理［12］对OMP算法的匹配过程和残差更新过程进行改进，提出基于粒子群优化的快速稀疏分解算法（fastPSO_OMP）。

2 基本理论

2.1 OMP算法的稀疏分解过程

利用稀疏分解对原始信号进行稀疏表示，能够得到信号的简洁表达，非正交分解能够对图像这类复杂信号进行有效稀疏表示，得到研究者们的关注。冗余字典的原子具有多样性，原子在信号组成的空间中足够的密。冗余字典标记为D={gγ}γ∈Γ，其中，构造原子的参数为γ，利用参数γ构造的原子定义为gγ，Γ为参数γ的集合。原始信号标记为f，维数为N。OMP算法实现信号稀疏分解的基本思想是：找到一些最优原子，这些最优原子的线性组合与原始信号的残差最小。具体步骤如下：

1）初始化残差r0=f，最优原子的索引集合Λ0=[]，OMP的分解次数k=1；

3）利用λk更新索引集合，Λk=Λk-1∪λk；

5）判断是否达到OMP的最大分解次数，若不满足则k=k+1，重复2）～4）；若满足则停止迭代。

经过K次迭代后，剩余残差的范数满足：

其中，gγk是第k次分解找到的最优原子。用最优原子表示原始信号的逼近，重构信号表示为

其中，ΘΛ表示最优原子构成的子矩阵。

2.2 粒子群优化和Hermitian求逆引理

1）粒子群优化

粒子群优化的目标是用粒子的进化运动代替冗余字典的原子遍历过程。假设搜索空间为L维，种群包含Pop个粒子，其中第i个粒子代表其在搜索空间中的位置，表示为向量Zi=(zi1，zi2，…，zil，…，ziL)，i=1，2，…，Pop，l=1，2，...，L。设定第i个粒子的 速 度 为， 个体极值为，种群的群体极值为。

粒子通过跟踪个体极值Pbest和群体极值Gbest更新其速度和位置［10］，即：

其中，iter为当前进化代数，w为惯性权重，c1和c2是非负常数，r1和r2是分布于［0，1］区间的随机数。

2）Hermitian求逆引理

Hermitian求逆引理能够实现Hermitian矩阵求逆的迭代求解，提高矩阵求逆速度。Hermitian矩阵的分块形式为

成立。

证明：令逆矩阵

由此导出：

由式（10）可得到：

将式（13）代入式（12）：

将式（14）代入式（13）：

将式（16）代入式（9），得到：

将式（14）～式（17）代入式（7），即得：

3 高光谱图像快速稀疏分解算法

3.1 Gabor字典

Gabor原子的生成函数［13］：

其 中 ，n=0，1，...，N，N为信号长度 ，是时频参数 ，其变化范围是：1≤s≤N，1≤u≤N，0≤υ≤2π，0≤ω≤2π。

3.2 算法实现过程

1）粒子群优化的实现过程

将OMP算法每次迭代中的匹配过程改进为粒子群优化对所有原子的寻优过程。粒子群优化中的粒子代表冗余字典的原子，通过对所有粒子的寻优，找到群体的极值位置，该极值位置对应的原子即为OMP每次迭代所寻找到的最优原子。

（1）粒子定义和适应度

粒子用四维向量来表示，定义为Zi=(s，u，υ，ω)，表征粒子好坏的参数是粒子代表的原子与残差的内积，考虑将粒子群算法的适应度函数表示为

（2）粒子群的初始化

（3）最优原子的寻优

2）Hermitian求逆引理的实现过程

与式（5）对比可知：

将Hk-1和Hk代入式（18），得到残差更新过程中矩阵求逆的递归表达式：

3.3 算法执行过程

本文提出的快速稀疏分解算法的流程图如图1所示，图中黑色粗实线框内是利用粒子群算法寻找最优原子的过程，黑色粗虚线框内是利用Hermitian求逆引理进行残差更新的过程。

算法的执行过程总结如下：

1）初始化残差r0=f，最优原子的集合Λ0=[]，fastPSO_OMP的分解次数k=1；

2）粒子群优化寻找最优原子：初始化粒子，计算粒子的适应度，确定Pbest和Gbest；根据速度和位置更新公式，粒子不断运动；更新Pbest和Gbest。检查是否达到粒子的最大进化代数Tmax，若达到则终止迭代，输出Gbest，否则继续运动；

3）利用Gbest更新索引集合，Λk=Λk-1∪Gbest；

5）判断是否满足迭代终止条件，若不满足则k=k+1，重复2）～4）；若满足则停止迭代。

利用PSO搜索得到的最优原子对原始信号进行稀疏表示，重构信号仍用式（2）计算。

图1 算法fastPSO_OMP的流程图

3.4 算法对比分析

算法fastPSO_OMP和OMP算法进行对比分析：

1）两种算法占用的计算机内存不同。OMP算法完成基于Gabor字典的稀疏分解时，需事先产生Gabor离散字典。假设信号长度N=256，原子个数将达到119756，占用计算机内存172MB。算法fastPSO_OMP则无需事先产生字典，粒子参数向量占用的内存基本可以忽略不计。

2）两种算法的稀疏表示能力不同。OMP算法中，Gabor字典事先确定，原子的时频参数完全固定，将捕捉不到特征复杂的图像信息。算法fastPSO_OMP的粒子能够在整个参数空间中运动，形成的原子比较灵活，能够生成更多类型的原子，适用于图像信号的稀疏分解，稀疏表示能力较固定的Gabor原子更强。

4 实验结果与分析

4.1 高光谱数据集

选择四组高光谱图像（数据来源：http：//aviris.jpl.nasa.gov）评估稀疏分解算法的性能。前两组数据是AVIRIS采集得到的场景Cuprite1和场景Cuprite2，数据集共有224个波段，去除异常波段和全零波段后，可用波段数为188，图像大小裁剪至256*256。第三组数据是由AVIRIS采集的场景Indian Pines，数据集共220个波段，去除水汽吸收波段后，可用波段数为200，图像大小裁剪至128*128。第四组数据是ROSIS采集的场景Pavia University，数据集包含115个波段，去除含噪波段后，可用波段数为103，图像大小裁剪至256*256。高光谱图像用二维矩阵X∈RNp×Nλ表示，其中，X=[X1，X2，…，Xb，…，XNλ]，Xb是第b个波段图像的向量化表示，Np是单个波段图像的像素点数，Nλ表示波段数。第50个波段的原始图像如图2所示。

4.2 实验参数设置

对图像进行分块处理，分块大小为B=16，则稀疏分解算法处理的信号长度为N=256。利用重构图像的峰值信噪比（Peak Signal-to-Noise Ratio，PSNR）和重构时间对算法性能进行评价。

单个图像块的重构PSNR的计算公式为

其中，Xb，bl是图像块的向量化表示，是重构图像块的向量化表示，是Xb，bl的峰值，是均方误差，定义为

图2 四组高光谱数据的第50个波段的原始图像

利用所提算法fastPSO_OMP进行稀疏分解，粒子群优化中最大进化代数Tmax的变化范围是5～50，间隔是5，种群大小Pop的变化范围是5～50，间隔是5，分解次数K的变化范围是10～100，间隔是10。图3给出的是Cuprite1重构图像的平均PSNR随参数的变化。当最大进化代数和种群个数逐渐增大时，粒子群优化需要从更多的粒子中经过多次进化才能找到最优原子，即在二者的影响下，时间复杂度均会不断增加。在保证重构精度的条件下，应尽量降低最大进化代数和种群数量。综合考虑，将所提算法fastPSO_OMP中粒子群优化的最大进化代数设置为Tmax=5，种群大小设置为Pop=10。

利用fastPSO_OMP算法和OMP算法对四组高光谱图像进行稀疏分解，设定算法在分解次数达到100次时终止。图4给出两种算法得到的重构PSNR随分解次数的变化，其中黑色粗线是OMP分解次数达到五十次时，利用50个最优原子对原始图像进行稀疏表示后，得到的重构图像精度。以Cuprite1和Indian Pines为例，将OMP算法利用50个原子得到的重构图像的PSNR作为标准，Cuprite1需要80个原子左右才能达到此标准，Indian Pines需要100个甚至更多的原子才能达到与OMP算法相同的精度。据此，将OMP算法设定为K=50，fastPSO_OMP算法设定为K=100。

图3 最大进化代数、种群大小和分解次数对重构的影响

图4 两种算法的重构PSNR随分解次数的变化

4.3 实验结果分析

利用OMP算法和fastPSO_OMP算法对四组高光谱图像进行稀疏分解，分析两种算法的性能。表1给出两种算法的实验结果，表中的平均PSNR是高光谱数据集所有波段图像重构PSNR的均值，运行时间是完成单个波段稀疏分解所需的时间，加速比是指相对比OMP算法的稀疏分解速度，所提算法fastPSO_OMP的加速倍数。

表1 OMP算法和fastPSO_OMP算法的性能对比

对于场景Cuprite1和Cuprite2，利用本文算法fastPSO_OMP得到的重构PSNR均与OMP算法得到的PSNR相当，而场景Indian Pines和Pavia University的重构PSNR略低于OMP算法的重构PSNR。实验结果能有效说明，所提算法能够从连续的时频参数范围内找到表示原始图像的最优原子，充分证明本文算法的可靠性。本文算法的执行效率远高于OMP算法，加速比能达到25倍左右。

5 结语

针对OMP算法计算效率低的问题，本文提出了基于粒子群优化的高光谱图像快速稀疏分解算法，采用粒子群思想优化OMP算法的匹配过程，并利用Hermitian求逆引理对其残差更新过程进行改进。四组高光谱图像的稀疏分解实验结果充分说明，所提算法fastPSO_OMP能够保证重构精度，并将OMP算法的计算效率提高约25倍。粒子群算法的随机性如何降低是值得进一步探究的问题。