王 麗,王 威,陳 博2(西安航空學(xué)院電子工程學(xué)院,西安70077)
2(咸陽師范學(xué)院物理與電子工程學(xué)院,陜西咸陽712000)
E-mail:wangli871016@163.com
高光譜圖像的空間和譜間信息較為豐富,其應(yīng)用范圍較廣[1-3].然而,高光譜圖像的數(shù)據(jù)量較大,給高光譜圖像帶來了挑戰(zhàn),增加了傳輸和存儲(chǔ)系統(tǒng)的壓力.因此,在研究高光譜圖像時(shí),尋找有效的壓縮技術(shù)對高光譜圖像的發(fā)展非常關(guān)鍵.
壓縮感知理論(Compressed Sensing,CS)[4,5]指出,若信號或圖像具有稀疏性或可壓縮性,可通過低于奈奎斯特速率進(jìn)行采樣,利用較少的測量值仍能精確獲取重構(gòu)信號.高光譜圖像的空間和譜間均存在冗余,具備可壓縮性[6].因此,利用壓縮感知理論可對高光譜圖像進(jìn)行處理,基于冗余字典和正交匹配追蹤算法(Orthogonal Matching Pursuit,OMP)[7]可實(shí)現(xiàn)高光譜圖像的采樣和重構(gòu).因高光譜圖像的波段數(shù)眾多,當(dāng)利用OMP算法實(shí)現(xiàn)基于冗余字典的重構(gòu)時(shí),其計(jì)算復(fù)雜度太高.
為提高OMP算法的計(jì)算效率,文獻(xiàn)[8]提出利用粒子群算法(Particle Swarm Optimization,PSO)[9]優(yōu)化 OMP 算法的匹配過程,實(shí)現(xiàn)高光譜圖像稀疏分解的最優(yōu)原子的搜索,但OMP算法的求逆過程仍然非常耗時(shí).因此本文不僅用PSO優(yōu)化OMP的匹配過程,同時(shí)利用Hermitian求逆引理[10]優(yōu)化殘差更新過程,實(shí)現(xiàn)殘差的遞歸更新,應(yīng)用于高光譜圖像的壓縮感知重構(gòu),提出改進(jìn)的粒子群優(yōu)化 OMP重構(gòu)算法(fastPSO_OMP).該算法能夠快速實(shí)現(xiàn)高光譜圖像的重構(gòu),計(jì)算復(fù)雜度低,且重構(gòu)精度與OMP算法相當(dāng).
假設(shè)x是長度為N的離散時(shí)域信號,其可在某組基Ψ=[ψ1,ψ2,…,ψN]上表示為:
其中,S是投影系數(shù)Sd=?x,ψd?構(gòu)成的N*1維的列向量.
如果S中大部分系數(shù)為零或者排序后呈指數(shù)級衰減并趨近于零,則認(rèn)為信號x是可壓縮的或在Ψ域上是稀疏的.用一個(gè)M*N維的觀測矩陣Φ(MN)對信號x進(jìn)行測量,得到M維的觀測值y:
其中,Θ=ΦΨ是一個(gè)M*N維的矩陣.壓縮感知理論指出,在已知y和Θ的情況下,可利用重構(gòu)算法恢復(fù)出稀疏信號S,最后由公式(1)得到原始信號x.
傳統(tǒng)OMP算法的重構(gòu)過程如下:
1)初始化:殘差r0=y,最優(yōu)原子的索引集合Λ0=[],迭代次數(shù)k=1;
2)遍歷矩陣Θ中的所有原子(列向量),利用內(nèi)積選取最佳匹配原子的索引,其中,Θj是矩陣Θ的第j列;
3)利用 λk更新索引集合,
5)判斷是否達(dá)到最大迭代步數(shù)K,若不滿足則k=k+1,重復(fù)2)~4);若滿足則停止迭代.
重構(gòu)信號為:
其中,ΨΛK表示稀疏基Ψ中對應(yīng)于索引集合ΛK的列向量構(gòu)成的子矩陣.
由上述重構(gòu)過程可知,OMP的匹配過程在每次迭代時(shí)均需要遍歷所有原子以尋找最佳匹配原子.因冗余字典原子個(gè)數(shù)多,其匹配計(jì)算量大;同時(shí)其殘差更新中存在大型矩陣的求逆運(yùn)算,非常耗時(shí).針對以上耗時(shí)過程,本文分別利用粒子群算法和Hermitian求逆引理對其進(jìn)行優(yōu)化.
粒子群算法的基本思想是通過群體中個(gè)體之間的協(xié)作和群體中信息的共享產(chǎn)生進(jìn)化優(yōu)勢,從而達(dá)到尋找最優(yōu)解的目的.本文利用粒子代表冗余字典的原子,用位置、速度和適應(yīng)度值表示.假設(shè)搜索空間為L維,種群Z=(Z1,Z2,…,Zm)包含m個(gè)粒子,第i個(gè)粒子表示為向量 Zi=(zi1,zi2,…,ziL),代表其在搜索空間中的位置.根據(jù)適應(yīng)度函數(shù)可計(jì)算出粒子Zi的適應(yīng)度值.第 i個(gè)粒子的速度為 Vi=(vi1,vi2,…,viL),個(gè)體極值為 Pbesti=(Pi1,Pi2,…,PiL),群體極值為 Gbest=(G1,G2,…,GL).粒子通過跟蹤個(gè)體極值Pbest和群體極值Gbest更新其速度和位置[9],隨著迭代的進(jìn)行,粒子不斷運(yùn)動(dòng)以找到群體的極值,即最優(yōu)原子,大大縮短運(yùn)算時(shí)間.
隨著迭代次數(shù)的增加,索引集字典的維度越來越高,導(dǎo)致求逆運(yùn)算-1的計(jì)算復(fù)雜度越來越高.因是正定Hermitian矩陣,本文利用Hermitian求逆引理實(shí)現(xiàn)矩陣的遞歸求逆,提高矩陣求逆速度.Hermitian求逆引理表示如下:
Hermitian矩陣的分塊形式為:
成立.即利用前次迭代的求逆結(jié)果線性表示當(dāng)次迭代的求逆結(jié)果,大大提高運(yùn)行效率.
本文利用Gabor冗余字典作為高光譜圖像的稀疏基,Gabor字典的生成函數(shù)表示為[11]:冗余度很高,原子個(gè)數(shù)為52(N log2N+N-1).
由于OMP算法的匹配過程和殘差更新過程均非常耗時(shí),本文提出一種改進(jìn)的粒子群優(yōu)化OMP重構(gòu)算法(fastPSO_OMP),利用粒子群算法的局部尋優(yōu)能力優(yōu)化OMP的匹配過程,快速找到最優(yōu)原子;并且利用Hermitian求逆引理加速OMP的殘差更新過程,實(shí)現(xiàn)殘差的遞歸更新,以降低算法的計(jì)算復(fù)雜度.下面分別介紹粒子群優(yōu)化和Hermitian求逆引理的實(shí)現(xiàn)過程.
3.2.1 粒子群優(yōu)化的實(shí)現(xiàn)過程
1)粒子位置和適應(yīng)度
粒子位置為(s,u,υ,ω),是一個(gè) 4維向量,選擇 OMP 的匹配函數(shù)作為適應(yīng)度函數(shù):
2)粒子群的初始化
假定粒子群的群體規(guī)模為m,初始群體為Z0=(Z01,Z02,…),第i個(gè)粒子的初始位置為
其中,n=0,1,…,N,N 為信號長度,γ =(s,u,υ,ω)是時(shí)頻參數(shù).時(shí)頻參數(shù)按以下方法離散化:γ =(aα1,α2aα1Δu,α3a-α1Δυ,α4Δω),其中,a=2,,其取值分別是(s,u,υ,ω)取值范圍中的一個(gè)隨機(jī)數(shù).根據(jù)生成函數(shù)計(jì)算粒子所生成的原子,并計(jì)算適應(yīng)度.第i個(gè)粒子的最優(yōu)位置為Pbesti=Z0i,初始群體極值為1,2,…,m.同時(shí)設(shè)置每個(gè)粒子的初始速度.
3)最優(yōu)原子的確定
在迭代過程中,更新每個(gè)粒子的速度和位置,并評價(jià)其適應(yīng)度.若 fitness()>Pbesti,則更新個(gè)體極值位置 Pbesti=,否則保持不變.若fitness()>Gbest,則更新群體極值位置Gbest=Ziiter,否則保持不變.迭代得到的群體極值位置所對應(yīng)的原子即為最優(yōu)原子.
3.2.2 Hermitian求逆引理的實(shí)現(xiàn)過程
fastPSO_OMP的每步迭代中,原子字典的更新表示為ΨΛk=[ΨΛk-1gγk],其中,ΨΛk-1是第 k -1 次迭代后形成的原子字典,gγk是第 k 次迭代找到的最優(yōu)原子,則 ΘΛk=[ΘΛk-1與公式(4)對比可知:
將Hk-1和Hk代入公式(5),得到殘差更新過程中矩陣求逆的遞歸表達(dá)式,見公式(11).通過遞歸方式實(shí)現(xiàn)傳統(tǒng)OMP算法的殘差更新過程,提高算法的執(zhí)行效率.
本文提出的改進(jìn)的粒子群優(yōu)化OMP重構(gòu)算法的流程圖如圖1所示,圖中黑色粗實(shí)線框內(nèi)是利用粒子群算法尋找最優(yōu)原子的過程,黑色粗虛線框內(nèi)是利用Hermitian求逆引理進(jìn)行殘差更新的過程.成的最優(yōu)原子gγk,根據(jù)公式(11),利用第k-1步迭代所得到的(-1遞歸求解-1,并計(jì)算殘差 rk=y
圖1 改進(jìn)的fastPSO_OMP重構(gòu)算法流程圖Fig.1 Framework of fastPSO_OMPreconstruction algorithm
5)判斷是否達(dá)到最大迭代次數(shù)K,若不滿足則k=k+1,重復(fù)2)~4);若滿足則停止迭代.
重構(gòu)信號表達(dá)式仍為公式(3),其中,ΨΛK表示第K步迭代后形成的原子字典.
利用OMP算法、PSO_OMP算法和fastPSO_OMP算法對四組高光譜圖像1數(shù)據(jù)來源:http://aviris.jpl.nasa.gov進(jìn)行壓縮感知重構(gòu).前兩組數(shù)據(jù)是AVIRIS采集得到的場景Cuprite1和場景Cuprite2,數(shù)據(jù)集共有224個(gè)波段,去除異常波段和全零波段后,可用波段數(shù)為188,圖像大小裁剪至256×256.第三組數(shù)據(jù)是由AVIRIS采集的場景Indian Pines,數(shù)據(jù)集共220個(gè)波段,去除水汽吸收波段后,可用波段數(shù)為200,圖像大小裁剪至128×128.第四組數(shù)據(jù)是ROSIS采集的場景Pavia University,數(shù)據(jù)集包含115個(gè)波段,去除含噪波段后,可用波段數(shù)為103,圖像大小裁剪至256×256.為保證數(shù)據(jù)的統(tǒng)一性,對高光譜圖像進(jìn)行處理時(shí),對原始采集得到的圖像均重新進(jìn)行量化,將每個(gè)波段圖像的量化位數(shù)設(shè)置為8位,即圖像像素值范圍是0~255.四組高光譜數(shù)據(jù)第50個(gè)波段的原始圖像如圖2所示.
因計(jì)算機(jī)內(nèi)存有限,實(shí)驗(yàn)中利用高斯隨機(jī)矩陣對高光譜圖像進(jìn)行分塊壓縮感知采樣,分塊大小為16,并利用OMP算法、PSO_OMP算法和fastPSO_OMP算法進(jìn)行重構(gòu).三種重構(gòu)算法的迭代次數(shù)K=50,粒子群算法的迭代步數(shù)Tmax=5,種群大小m=10.
三種重構(gòu)算法重構(gòu)單個(gè)波段所需的匹配時(shí)間和更新時(shí)間見表1.從表中可以看出,僅用粒子群優(yōu)化OMP匹配過程的PSO_OMP算法,可以將匹配時(shí)間降低兩個(gè)數(shù)量級,但對更新時(shí)間的影響較小.本文進(jìn)一步用Hermitian求逆引理對OMP的更新過程進(jìn)行改進(jìn),提出的fastPSO_OMP算法能夠?qū)MP算法的更新時(shí)間降低一個(gè)數(shù)量級.加速原因在于:1)OMP需要遍歷Gabor字典中的每個(gè)原子,而PSO僅需搜索少量的空間參數(shù)點(diǎn)就能尋找到最優(yōu)原子;具體來說,在每次迭代時(shí),OMP算法需遍歷119756個(gè)原子才能找到一個(gè)最優(yōu)原子,而PSO算法僅需搜索50個(gè)空間參數(shù)點(diǎn)(迭代次數(shù)*種群數(shù)量)就能找到一個(gè)最優(yōu)原子,計(jì)算量大大降低.2)引入Hermitian求逆引理,每次迭代過程中僅需少量矩陣計(jì)算,采用遞歸方式實(shí)現(xiàn)矩陣求逆,運(yùn)算量大大降低,算法運(yùn)行效率提高.
表1 三種算法重構(gòu)單個(gè)波段所需的匹配時(shí)間和更新時(shí)間Table 1 Matching time and update time of three reconstruction algorithms for reconstructing one band
表2 三種算法的平均PSNR及加速比Table 2 Average PSNR and accelerating times of three reconstruction algorithms for reconstructing one band
利用重構(gòu)圖像的峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)、運(yùn)行時(shí)間(runtime)及重構(gòu)圖像的主觀視覺對算法性能進(jìn)行評價(jià).不同采樣率下,所有波段的平均PSNR及算法的加速比見表2.表中的平均PSNR是將每個(gè)波段的PSNR對所有波段進(jìn)行平均得到的.表中的加速比是指將OMP算法的計(jì)算速度設(shè)為1,PSO_OMP算法和本文算法fastPSO_OMP的計(jì)算速度相對于OMP算法速度的倍數(shù).從表中可以看出,對于場景Cuprite1和場景Cuprite2來說,算法PSO_OMP和算法fastPSO_OMP的PSNR與OMP算法相當(dāng),在采樣率為0.1和0.5時(shí)有提高的優(yōu)勢,其他采樣率下的差值能保證在1dB以內(nèi).對于場景Indian Pines和場景Pavia University來說,算法PSO_OMP和所提算法fastPSO_OMP的重構(gòu)性能在五個(gè)采樣率下均有不同程度的提高,特別是在采樣率為0.1時(shí),提高值分別為1.76dB和1.31dB,算法性能得到保證.同時(shí)從表中可以看出,相較于算法PSO_OMP,本文算法fastPSO_OMP進(jìn)一步地提高了執(zhí)行效率,加速比最高可達(dá)到18.8657.
重構(gòu)性能提高的原因在于:利用PSO算法搜索最優(yōu)原子時(shí),粒子位置在整個(gè)參數(shù)空間中變化,能夠搜索得到最能表示圖像特征的原子,而OMP算法無法在事先構(gòu)造的冗余字典中找到與信號特征最相符的原子,因此利用粒子群算法能夠在一定程度上提高重構(gòu)圖像的質(zhì)量.執(zhí)行效率提高的原因在于:利用PSO搜索最優(yōu)原子時(shí),僅需搜索50個(gè)參數(shù)空間點(diǎn),其運(yùn)行時(shí)間遠(yuǎn)遠(yuǎn)低于OMP算法的匹配過程;同時(shí)本文引入Hermitian求逆引理,采用遞歸方式解決了殘差更新的矩陣求逆問題,降低了OMP算法中大型矩陣的求逆時(shí)間.故本文算法能夠在較少時(shí)間內(nèi)完成高光譜圖像壓縮感知重構(gòu),在保證 重構(gòu)精度的條件下,大大提高計(jì)算效率.個(gè)波段,采樣率為0.5.算法PSO_OMP和本文算法fastPSO_OMP均能夠得到重構(gòu)質(zhì)量優(yōu)于OMP算法的圖像.
圖2 4組高光譜數(shù)據(jù)的第50個(gè)波段的原始圖像Fig.2 50th band original images of four hyperspectral datasets
圖3 數(shù)據(jù)集Cuprite2的重構(gòu)圖像與原始圖像對比Fig.3 Comparison of reconstructed image of dataset Cuprite2 with its original image
利用OMP算法、PSO_OMP算法和fastPSO_OMP算法得到的重構(gòu)Cuprite2圖像如圖3所示,圖中為重構(gòu)圖像的第40
為提高OMP算法處理基于冗余字典的高光譜圖像壓縮感知重構(gòu)問題的計(jì)算效率,本文采用粒子群算法優(yōu)化OMP的匹配過程,并引入Hermitian求逆引理加速其殘差更新過程,提出改進(jìn)的粒子群優(yōu)化OMP重構(gòu)算法(fastPSO_OMP).將所提算法應(yīng)用于高光譜圖像壓縮感知重構(gòu),實(shí)驗(yàn)結(jié)果表明,與OMP算法相比,本文算法能夠在保證重構(gòu)精度的條件下,大大提高計(jì)算效率.