王 麗 孫長(zhǎng)杰王 威

（1.西安航空學(xué)院 西安 710077）（2.中國(guó)飛行試驗(yàn)研究院 西安 710089）

1 引言

高光譜圖像（Hyperspectral Images，HSIs）包含豐富的空間幾何信息和光譜特征信息，適用于戰(zhàn)場(chǎng)地形偵察［1］、目標(biāo)探測(cè)與識(shí)別［2～3］、分類［4～5］等軍事和民用領(lǐng)域。高光譜成像光譜儀的研制［6］、高光譜圖像解譯研究［7］，已成為備受關(guān)注的前沿?zé)狳c(diǎn)。應(yīng)用領(lǐng)域的不斷擴(kuò)展，要求高光譜圖像能夠提供更細(xì)致的目標(biāo)信息。問(wèn)題是，空間和光譜分辨率越高，對(duì)系統(tǒng)的數(shù)據(jù)存儲(chǔ)和傳輸能力的挑戰(zhàn)更高。將信號(hào)稀疏分解理論用于高光譜圖像［8］，便于對(duì)高光譜圖像進(jìn)行后續(xù)處理。高光譜圖像波段多，如何尋找計(jì)算效率高的稀疏分解算法成為研究的關(guān)鍵。

經(jīng)典的稀疏分解算法是正交匹配追蹤（Orthog-onal Matching Pursuit，OMP）［9］算法，該算法能夠從原子庫(kù)中依次查找到能夠?qū)π盘?hào)進(jìn)行稀疏表示的原子，重構(gòu)精度高。高光譜圖像特征復(fù)雜，采用冗余度高的冗余字典才能完成其稀疏表示。正交匹配追蹤算法在處理此類問(wèn)題時(shí)，計(jì)算時(shí)間在現(xiàn)有計(jì)算條件下無(wú)法忍受。鑒于進(jìn)化算法［10］尋優(yōu)能力強(qiáng)的特點(diǎn)，本文采用粒子群算法（Particle Swarm Optimization，PSO）［11］和Hermitian求逆引理［12］對(duì)OMP算法的匹配過(guò)程和殘差更新過(guò)程進(jìn)行改進(jìn)，提出基于粒子群優(yōu)化的快速稀疏分解算法（fastPSO_OMP）。

2 基本理論

2.1 OMP算法的稀疏分解過(guò)程

利用稀疏分解對(duì)原始信號(hào)進(jìn)行稀疏表示，能夠得到信號(hào)的簡(jiǎn)潔表達(dá)，非正交分解能夠?qū)D像這類復(fù)雜信號(hào)進(jìn)行有效稀疏表示，得到研究者們的關(guān)注。冗余字典的原子具有多樣性，原子在信號(hào)組成的空間中足夠的密。冗余字典標(biāo)記為D={gγ}γ∈Γ，其中，構(gòu)造原子的參數(shù)為γ，利用參數(shù)γ構(gòu)造的原子定義為gγ，Γ為參數(shù)γ的集合。原始信號(hào)標(biāo)記為f，維數(shù)為N。OMP算法實(shí)現(xiàn)信號(hào)稀疏分解的基本思想是：找到一些最優(yōu)原子，這些最優(yōu)原子的線性組合與原始信號(hào)的殘差最小。具體步驟如下：

1）初始化殘差r0=f，最優(yōu)原子的索引集合Λ0=[]，OMP的分解次數(shù)k=1；

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

5）判斷是否達(dá)到OMP的最大分解次數(shù)，若不滿足則k=k+1，重復(fù)2）～4）；若滿足則停止迭代。

經(jīng)過(guò)K次迭代后，剩余殘差的范數(shù)滿足：

其中，gγk是第k次分解找到的最優(yōu)原子。用最優(yōu)原子表示原始信號(hào)的逼近，重構(gòu)信號(hào)表示為

其中，ΘΛ表示最優(yōu)原子構(gòu)成的子矩陣。

2.2 粒子群優(yōu)化和Hermitian求逆引理

1）粒子群優(yōu)化

粒子群優(yōu)化的目標(biāo)是用粒子的進(jìn)化運(yùn)動(dòng)代替冗余字典的原子遍歷過(guò)程。假設(shè)搜索空間為L(zhǎng)維，種群包含Pop個(gè)粒子，其中第i個(gè)粒子代表其在搜索空間中的位置，表示為向量Zi=(zi1，zi2，…，zil，…，ziL)，i=1，2，…，Pop，l=1，2，...，L。設(shè)定第i個(gè)粒子的 速 度 為， 個(gè)體極值為，種群的群體極值為。

粒子通過(guò)跟蹤個(gè)體極值Pbest和群體極值Gbest更新其速度和位置［10］，即：

其中，iter為當(dāng)前進(jìn)化代數(shù)，w為慣性權(quán)重，c1和c2是非負(fù)常數(shù)，r1和r2是分布于［0，1］區(qū)間的隨機(jī)數(shù)。

2）Hermitian求逆引理

Hermitian求逆引理能夠?qū)崿F(xiàn)Hermitian矩陣求逆的迭代求解，提高矩陣求逆速度。Hermitian矩陣的分塊形式為

成立。

證明：令逆矩陣

由此導(dǎo)出：

由式（10）可得到：

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

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

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

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

3 高光譜圖像快速稀疏分解算法

3.1 Gabor字典

Gabor原子的生成函數(shù)［13］：

其 中 ，n=0，1，...，N，N為信號(hào)長(zhǎng)度 ，是時(shí)頻參數(shù) ，其變化范圍是：1≤s≤N，1≤u≤N，0≤υ≤2π，0≤ω≤2π。

3.2 算法實(shí)現(xiàn)過(guò)程

1）粒子群優(yōu)化的實(shí)現(xiàn)過(guò)程

將OMP算法每次迭代中的匹配過(guò)程改進(jìn)為粒子群優(yōu)化對(duì)所有原子的尋優(yōu)過(guò)程。粒子群優(yōu)化中的粒子代表冗余字典的原子，通過(guò)對(duì)所有粒子的尋優(yōu)，找到群體的極值位置，該極值位置對(duì)應(yīng)的原子即為OMP每次迭代所尋找到的最優(yōu)原子。

（1）粒子定義和適應(yīng)度

粒子用四維向量來(lái)表示，定義為Zi=(s，u，υ，ω)，表征粒子好壞的參數(shù)是粒子代表的原子與殘差的內(nèi)積，考慮將粒子群算法的適應(yīng)度函數(shù)表示為

（2）粒子群的初始化

（3）最優(yōu)原子的尋優(yōu)

2）Hermitian求逆引理的實(shí)現(xiàn)過(guò)程

與式（5）對(duì)比可知：

將Hk-1和Hk代入式（18），得到殘差更新過(guò)程中矩陣求逆的遞歸表達(dá)式：

3.3 算法執(zhí)行過(guò)程

本文提出的快速稀疏分解算法的流程圖如圖1所示，圖中黑色粗實(shí)線框內(nèi)是利用粒子群算法尋找最優(yōu)原子的過(guò)程，黑色粗虛線框內(nèi)是利用Hermitian求逆引理進(jìn)行殘差更新的過(guò)程。

算法的執(zhí)行過(guò)程總結(jié)如下：

1）初始化殘差r0=f，最優(yōu)原子的集合Λ0=[]，fastPSO_OMP的分解次數(shù)k=1；

2）粒子群優(yōu)化尋找最優(yōu)原子：初始化粒子，計(jì)算粒子的適應(yīng)度，確定Pbest和Gbest；根據(jù)速度和位置更新公式，粒子不斷運(yùn)動(dòng)；更新Pbest和Gbest。檢查是否達(dá)到粒子的最大進(jìn)化代數(shù)Tmax，若達(dá)到則終止迭代，輸出Gbest，否則繼續(xù)運(yùn)動(dòng)；

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

5）判斷是否滿足迭代終止條件，若不滿足則k=k+1，重復(fù)2）～4）；若滿足則停止迭代。

利用PSO搜索得到的最優(yōu)原子對(duì)原始信號(hào)進(jìn)行稀疏表示，重構(gòu)信號(hào)仍用式（2）計(jì)算。

圖1 算法fastPSO_OMP的流程圖

3.4 算法對(duì)比分析

算法fastPSO_OMP和OMP算法進(jìn)行對(duì)比分析：

1）兩種算法占用的計(jì)算機(jī)內(nèi)存不同。OMP算法完成基于Gabor字典的稀疏分解時(shí)，需事先產(chǎn)生Gabor離散字典。假設(shè)信號(hào)長(zhǎng)度N=256，原子個(gè)數(shù)將達(dá)到119756，占用計(jì)算機(jī)內(nèi)存172MB。算法fastPSO_OMP則無(wú)需事先產(chǎn)生字典，粒子參數(shù)向量占用的內(nèi)存基本可以忽略不計(jì)。

2）兩種算法的稀疏表示能力不同。OMP算法中，Gabor字典事先確定，原子的時(shí)頻參數(shù)完全固定，將捕捉不到特征復(fù)雜的圖像信息。算法fastPSO_OMP的粒子能夠在整個(gè)參數(shù)空間中運(yùn)動(dòng)，形成的原子比較靈活，能夠生成更多類型的原子，適用于圖像信號(hào)的稀疏分解，稀疏表示能力較固定的Gabor原子更強(qiáng)。

4 實(shí)驗(yàn)結(jié)果與分析

4.1 高光譜數(shù)據(jù)集

選擇四組高光譜圖像（數(shù)據(jù)來(lái)源：http：//aviris.jpl.nasa.gov）評(píng)估稀疏分解算法的性能。前兩組數(shù)據(jù)是AVIRIS采集得到的場(chǎng)景Cuprite1和場(chǎng)景Cuprite2，數(shù)據(jù)集共有224個(gè)波段，去除異常波段和全零波段后，可用波段數(shù)為188，圖像大小裁剪至256*256。第三組數(shù)據(jù)是由AVIRIS采集的場(chǎng)景Indian Pines，數(shù)據(jù)集共220個(gè)波段，去除水汽吸收波段后，可用波段數(shù)為200，圖像大小裁剪至128*128。第四組數(shù)據(jù)是ROSIS采集的場(chǎng)景Pavia University，數(shù)據(jù)集包含115個(gè)波段，去除含噪波段后，可用波段數(shù)為103，圖像大小裁剪至256*256。高光譜圖像用二維矩陣X∈RNp×Nλ表示，其中，X=[X1，X2，…，Xb，…，XNλ]，Xb是第b個(gè)波段圖像的向量化表示，Np是單個(gè)波段圖像的像素點(diǎn)數(shù)，Nλ表示波段數(shù)。第50個(gè)波段的原始圖像如圖2所示。

4.2 實(shí)驗(yàn)參數(shù)設(shè)置

對(duì)圖像進(jìn)行分塊處理，分塊大小為B=16，則稀疏分解算法處理的信號(hào)長(zhǎng)度為N=256。利用重構(gòu)圖像的峰值信噪比（Peak Signal-to-Noise Ratio，PSNR）和重構(gòu)時(shí)間對(duì)算法性能進(jìn)行評(píng)價(jià)。

單個(gè)圖像塊的重構(gòu)PSNR的計(jì)算公式為

其中，Xb，bl是圖像塊的向量化表示，是重構(gòu)圖像塊的向量化表示，是Xb，bl的峰值，是均方誤差，定義為

圖2 四組高光譜數(shù)據(jù)的第50個(gè)波段的原始圖像

利用所提算法fastPSO_OMP進(jìn)行稀疏分解，粒子群優(yōu)化中最大進(jìn)化代數(shù)Tmax的變化范圍是5～50，間隔是5，種群大小Pop的變化范圍是5～50，間隔是5，分解次數(shù)K的變化范圍是10～100，間隔是10。圖3給出的是Cuprite1重構(gòu)圖像的平均PSNR隨參數(shù)的變化。當(dāng)最大進(jìn)化代數(shù)和種群個(gè)數(shù)逐漸增大時(shí)，粒子群優(yōu)化需要從更多的粒子中經(jīng)過(guò)多次進(jìn)化才能找到最優(yōu)原子，即在二者的影響下，時(shí)間復(fù)雜度均會(huì)不斷增加。在保證重構(gòu)精度的條件下，應(yīng)盡量降低最大進(jìn)化代數(shù)和種群數(shù)量。綜合考慮，將所提算法fastPSO_OMP中粒子群優(yōu)化的最大進(jìn)化代數(shù)設(shè)置為T(mén)max=5，種群大小設(shè)置為Pop=10。

利用fastPSO_OMP算法和OMP算法對(duì)四組高光譜圖像進(jìn)行稀疏分解，設(shè)定算法在分解次數(shù)達(dá)到100次時(shí)終止。圖4給出兩種算法得到的重構(gòu)PSNR隨分解次數(shù)的變化，其中黑色粗線是OMP分解次數(shù)達(dá)到五十次時(shí)，利用50個(gè)最優(yōu)原子對(duì)原始圖像進(jìn)行稀疏表示后，得到的重構(gòu)圖像精度。以Cuprite1和Indian Pines為例，將OMP算法利用50個(gè)原子得到的重構(gòu)圖像的PSNR作為標(biāo)準(zhǔn)，Cuprite1需要80個(gè)原子左右才能達(dá)到此標(biāo)準(zhǔn)，Indian Pines需要100個(gè)甚至更多的原子才能達(dá)到與OMP算法相同的精度。據(jù)此，將OMP算法設(shè)定為K=50，fastPSO_OMP算法設(shè)定為K=100。

圖3 最大進(jìn)化代數(shù)、種群大小和分解次數(shù)對(duì)重構(gòu)的影響

圖4 兩種算法的重構(gòu)PSNR隨分解次數(shù)的變化

4.3 實(shí)驗(yàn)結(jié)果分析

利用OMP算法和fastPSO_OMP算法對(duì)四組高光譜圖像進(jìn)行稀疏分解，分析兩種算法的性能。表1給出兩種算法的實(shí)驗(yàn)結(jié)果，表中的平均PSNR是高光譜數(shù)據(jù)集所有波段圖像重構(gòu)PSNR的均值，運(yùn)行時(shí)間是完成單個(gè)波段稀疏分解所需的時(shí)間，加速比是指相對(duì)比OMP算法的稀疏分解速度，所提算法fastPSO_OMP的加速倍數(shù)。

表1 OMP算法和fastPSO_OMP算法的性能對(duì)比

對(duì)于場(chǎng)景Cuprite1和Cuprite2，利用本文算法fastPSO_OMP得到的重構(gòu)PSNR均與OMP算法得到的PSNR相當(dāng)，而場(chǎng)景Indian Pines和Pavia University的重構(gòu)PSNR略低于OMP算法的重構(gòu)PSNR。實(shí)驗(yàn)結(jié)果能有效說(shuō)明，所提算法能夠從連續(xù)的時(shí)頻參數(shù)范圍內(nèi)找到表示原始圖像的最優(yōu)原子，充分證明本文算法的可靠性。本文算法的執(zhí)行效率遠(yuǎn)高于OMP算法，加速比能達(dá)到25倍左右。

5 結(jié)語(yǔ)

針對(duì)OMP算法計(jì)算效率低的問(wèn)題，本文提出了基于粒子群優(yōu)化的高光譜圖像快速稀疏分解算法，采用粒子群思想優(yōu)化OMP算法的匹配過(guò)程，并利用Hermitian求逆引理對(duì)其殘差更新過(guò)程進(jìn)行改進(jìn)。四組高光譜圖像的稀疏分解實(shí)驗(yàn)結(jié)果充分說(shuō)明，所提算法fastPSO_OMP能夠保證重構(gòu)精度，并將OMP算法的計(jì)算效率提高約25倍。粒子群算法的隨機(jī)性如何降低是值得進(jìn)一步探究的問(wèn)題。