李 卓，徐 哲，陳 昕，李淑琴

(1.網(wǎng)絡(luò)文化與數(shù)字傳播北京市重點(diǎn)實(shí)驗(yàn)室(北京信息科技大學(xué))，北京 100101； 2.北京信息科技大學(xué) 計(jì)算機(jī)學(xué)院，北京 100101)(*通信作者電子郵箱lizhuo@bistu.edu.cn)

0 引言

高光譜圖像具有高空間分辨率和高光譜分辨率以及較多的光譜通道數(shù)，可以從空間對(duì)地觀測(cè)的同時(shí)獲取眾多連續(xù)波段的地物光譜圖像，達(dá)到從空間直接識(shí)別地球表面物質(zhì)、識(shí)別偽裝的目的[1]。高光譜圖像光譜范圍覆蓋從可見(jiàn)光到紅外，波段數(shù)量達(dá)到數(shù)百個(gè)波段[2]，光譜分辨率為10～l00 nm，是一種三維立體圖像，即在普通二維圖像的基礎(chǔ)上又多了一維光譜信息，因此高光譜遙感圖像的數(shù)據(jù)量龐大。高光譜成像儀器產(chǎn)生的實(shí)時(shí)數(shù)據(jù)量可達(dá)到每秒幾百M(fèi)b。

在軍事、醫(yī)療等領(lǐng)域，光譜成像技術(shù)已得到了越來(lái)越廣泛的應(yīng)用。如果在成像端不對(duì)原始數(shù)據(jù)進(jìn)行壓縮處理，不可能滿足實(shí)時(shí)傳輸場(chǎng)景的要求，因此，對(duì)高光譜遙感圖像壓縮編碼的研究具有重要的應(yīng)用價(jià)值。

與普通圖像相比，高光譜圖像的主要特點(diǎn)有：1)圖像的空間相關(guān)性較低，譜間相關(guān)性較強(qiáng)。由于高光譜圖像的空間結(jié)構(gòu)特征很難有整體上的一致性，故像素值連續(xù)性較差，空間相關(guān)性較低。譜間相關(guān)性是指在同一空間位置上，各個(gè)譜段的光譜圖像像素具有相似性。產(chǎn)生這種相似性的主要原因是同一目標(biāo)在連續(xù)的波段上的反射值也具有連續(xù)性。較強(qiáng)的譜間相關(guān)性是高光譜圖像壓縮的主要依據(jù)。2)壓縮保真度要求高。由于需要從壓縮后的圖像序列中提取光譜信息，解壓縮后的圖像保真度與光譜分析準(zhǔn)確度成正相關(guān)，因而對(duì)圖像壓縮方法提出了很高的要求。

高光譜圖像的壓縮類型主要分為無(wú)損和有損兩類。無(wú)損壓縮主要應(yīng)用在對(duì)數(shù)據(jù)保真度要求極高的場(chǎng)景下，往往壓縮比較低，在3∶1至5∶1之間。2012年，國(guó)際空間數(shù)據(jù)系統(tǒng)咨詢委員會(huì)(Consultative Committee for Space Data Systems， CCSDS)發(fā)布了高光譜圖像數(shù)據(jù)無(wú)損壓縮標(biāo)準(zhǔn)CCSDS 123.0-B-1[3]。有損壓縮在降低了數(shù)據(jù)保真度的同時(shí)，壓縮比也可實(shí)現(xiàn)大幅提高，可達(dá)到64∶1。

高光譜圖像的有損壓縮算法可以分為三類：1)將經(jīng)典的二維壓縮算法擴(kuò)展到三維空間，如3D-SPIHT(3-D Set Partitioning In Hierarchical Trees)算法[4]和3D-SPECK(3-D Set Partitioning Embedded bloCK)算法[5]等。2)基于靜態(tài)圖像壓縮標(biāo)準(zhǔn)的高光譜圖像壓縮算法，該類算法通常先采用主成分分析(Principal Component Analysis，PCA)或離散小波變換(Discrete Wavelet Transform， DWT)特征提取技術(shù)對(duì)高光譜圖像進(jìn)行譜間向量去相關(guān)，再利用靜態(tài)圖像壓縮標(biāo)準(zhǔn)對(duì)其空間向量進(jìn)行編碼[6]。在靜態(tài)圖像壓縮領(lǐng)域，對(duì)去空間冗余和量化編碼等技術(shù)已有了比較多的研究，比如JPEG2000[7]、WebP[8]、BPG(Better Portable Graphics)[9]等。JPEG2000基于小波變換和位平面編碼；WebP使用VP8視頻壓縮標(biāo)準(zhǔn)中的I幀編碼算法；BPG使用了HEVC的I幀編碼算法。通常情況下，PCA譜間去相關(guān)與JPEG2000相結(jié)合的壓縮算法(PCA+JPEG2000)比DWT譜間去相關(guān)與JPEG2000相結(jié)合的壓縮方法(DWT+JPEG2000)具有更優(yōu)的性能[6]。文獻(xiàn)[10]將視頻編碼標(biāo)準(zhǔn)H.264/AVC應(yīng)用于高光譜圖像，證明了基于H.264/AVC方法設(shè)計(jì)低功耗、實(shí)時(shí)的高光譜數(shù)據(jù)編碼器的可能性。3)基于張量的壓縮算法。張量可以理解為多維空間中的數(shù)組。傳統(tǒng)的圖像壓縮方法是基于矢量數(shù)據(jù)壓縮處理，改變了高維數(shù)據(jù)的空間結(jié)構(gòu)，而張量克服了這個(gè)缺陷，從而很好地保持了數(shù)據(jù)的空間結(jié)構(gòu)。Veganzones等[11]提出了一種面向高光譜數(shù)據(jù)壓縮的非負(fù)張量CP分解(Canonical Polyadic Decomposition)算法，旨在減少內(nèi)存占用和加速計(jì)算。高光譜圖像中含有大量的光譜變異性，因此，利用乘法更新算法生成張量的方式具有極高的計(jì)算復(fù)雜度，有時(shí)甚至不可行，目前在嵌入式實(shí)時(shí)壓縮環(huán)境還很難得到應(yīng)用。

去譜間冗余的方法主要有KLT(Karhunen-Loeve Transform)[12]、離散余弦變換(Discrete Cosine Transform，DCT)[13]、DWT[14]等。其中，KLT理論上是最佳的，但變換的基函數(shù)需要計(jì)算原始數(shù)據(jù)相關(guān)的協(xié)方差矩陣和特征矢量，導(dǎo)致計(jì)算量較大。以往受硬件計(jì)算性能的限制，KLT很難滿足實(shí)時(shí)壓縮場(chǎng)景的要求。隨著GPU在通用計(jì)算上的廣泛應(yīng)用，以及硬件計(jì)算性能的顯著發(fā)展，使得KLT實(shí)時(shí)壓縮場(chǎng)景中的應(yīng)用已成為可能。

本文基于KLT和HEVC設(shè)計(jì)了一種高光譜圖像實(shí)時(shí)壓縮算法，并利用NVIDIA Jetson TX1平臺(tái)設(shè)計(jì)并實(shí)現(xiàn)了CPU和GPU異構(gòu)并行壓縮處理系統(tǒng)。實(shí)驗(yàn)結(jié)果表明，該算法與DWT+JPEG2000[15]相比，在峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)方面提高了1.36 dB，延遲在100 ms以內(nèi)，能滿足實(shí)時(shí)性要求。

1 基于KLT和HEVC的高光譜圖像壓縮

基于KLT和HEVC的高光譜壓縮算法的整體工作流程見(jiàn)圖1。詳細(xì)的KLT和HEVC壓縮算法流程見(jiàn)圖2。

圖1 基于KLT和HEVC的TX1嵌入式

圖2 基于KLT和HEVC的高光譜實(shí)時(shí)壓縮算法流程

從圖1可以看出，主要的計(jì)算內(nèi)容包括KLT和HEVC幀內(nèi)壓縮，KLT負(fù)責(zé)去除譜間相關(guān)性，而HEVC負(fù)責(zé)去除空間相關(guān)性和編碼。

KLT已被證明是統(tǒng)計(jì)特征基礎(chǔ)上的最優(yōu)線性正交變換，經(jīng)過(guò)KLT后的系數(shù)之間互不相關(guān)，與DCT和DWT相比，數(shù)據(jù)能量更為集中，更有利于后續(xù)的壓縮處理。KLT與DCT以及DWT的不同之處在于其變換的基函數(shù)是不固定的，它依賴于待變換數(shù)據(jù)的統(tǒng)計(jì)特征。KLT需要對(duì)每組數(shù)據(jù)計(jì)算相應(yīng)的基函數(shù)，這也使得它對(duì)待變換數(shù)據(jù)具有最好的匹配效果，因此，KLT能夠獲得優(yōu)于其他變換的去相關(guān)性能。正因?yàn)镵LT針對(duì)每一組數(shù)據(jù)都需要重新計(jì)算它的基函數(shù)，導(dǎo)致KLT的計(jì)算量遠(yuǎn)遠(yuǎn)高于DCT和DWT。以往嵌入式平臺(tái)的計(jì)算能力無(wú)法支持KLT進(jìn)行實(shí)時(shí)計(jì)算的要求。2015年，NVIDIA發(fā)布了具有高并行處理能力的GPU嵌入式平臺(tái)Jetson TX1[16]，使得基于KLT的實(shí)時(shí)壓縮嵌入式平臺(tái)成為可能，因此本文選用了KLT作為去光譜相關(guān)的計(jì)算方法。

經(jīng)KLT去除譜間相關(guān)性后，可將變換后的數(shù)據(jù)看作一組互不相關(guān)的靜態(tài)圖像，使用靜態(tài)圖像壓縮算法去除譜間相關(guān)性。目前成熟的靜態(tài)圖像壓縮方法有JPEG、JPEG2000，針對(duì)視頻的編碼算法H.264和HEVC也支持靜態(tài)圖像壓縮。根據(jù)性能比較，HEVC的靜態(tài)圖像壓縮性能具有最優(yōu)性能，因此本文選用了HEVC作為去空間相關(guān)性和編碼部分的算法。

1.1 基于KLT的去譜間冗余

KLT去除譜間相關(guān)性主要思想是首先在圖像序列之間作一維變換：對(duì)變換后的前幾幅圖像，因?yàn)榧辛藞D像序列的主要能量，故進(jìn)行低倍數(shù)壓縮；而對(duì)后幾幅圖像，因?yàn)槠淠芰枯^低，故進(jìn)行高倍數(shù)壓縮。

KLT是從消除圖像間相關(guān)性均方差最小的意義上講的最佳變換，因此在理論上是去除譜間相關(guān)性的最優(yōu)方法。

首先介紹KLT消除譜間冗余的原理。假設(shè)原始數(shù)據(jù)高光譜圖像為：

(1)

令my表示向量y的平均值，y的協(xié)方差矩陣記為Cy,通過(guò)變換y=Kx,有：

my=E(y)=E(Kx)=KE(x)=Kmx

E[(Kx)(Kx)T]-(Kmx)(Kmx)T=

KCxKT=Λ=diag[λ0,λ1,…,λL-1]

(2)

寫成矩陣形式：

Cy=KCxKT=Λ=

(3)

從式(3)可以看出，光譜數(shù)據(jù)X經(jīng)過(guò)KLT后，其協(xié)方差矩陣為對(duì)角矩陣，即除了自身的方差之外，各波段之間的協(xié)方差均為零，這與波段之間的相關(guān)性為零等價(jià),因此KLT消除了光譜數(shù)據(jù)X的譜間相關(guān)性。

下面介紹如何構(gòu)造KLT矩陣：

1)計(jì)算協(xié)方差矩陣：

(4)

其中：i,j∈[1,B],Xk,i是向量Xk的第i個(gè)分量，mi是均值向量m的第i個(gè)分量。

2)計(jì)算矩陣C的特征向量ei(i=1,2,…,B)，得到變換矩陣K=(e1,e2,…,eB)T。

因此，向量X的KLT變換為Y=KX。由于K-1=KT,則相應(yīng)的逆變換為X=KTY。

在去譜間冗余方面，將原始高光譜圖像變換為二維矩陣，一維為光譜維，另一維為空間維。即將空間的二維轉(zhuǎn)換為一維形式。然后，將具有B個(gè)光譜波段的二維矩陣作為B維向量進(jìn)行KLT。變換后的新矩陣中的值越靠前，能量越高，即重要性越大。

1.2 面向GPU的KLT算法

高光譜圖像數(shù)據(jù)量較大，單幅圖像達(dá)到100 MB以上，每一幅高光譜圖像存在較多的譜段，譜段數(shù)量在一百個(gè)以上。由于單幅圖像數(shù)據(jù)量較大，運(yùn)用KLT對(duì)高光譜圖像去譜間冗余時(shí)，在CPU上處理速度較慢，因此考慮在GPU上實(shí)現(xiàn)KLT，使用arrayFire[17]的GPU上層封裝庫(kù)對(duì)圖像矩陣部分進(jìn)行運(yùn)算加速。算法的具體流程如圖3。

圖3 面向GPU的KLT算法流程

面向GPU的KLT算法偽代碼如算法1所示。

算法1 面向GPU的KLT算法。

輸入 高光譜圖像存放目錄路徑Dir，高光譜圖像名稱ImageName。

輸出 去譜間冗余后的高光譜圖像Y。

1)

Image ← read(Dir+ImageName)

2)

3)

ImageArray←ImageData-ImageMean

4)

ImageArrayT← transpose(ImageArray)

5)

CovArray←

matmul(ImageArrayT,ImageArray)/

(ImageArray.dims(0)-1)

//matmul為arrayfire內(nèi)置矩陣乘法函數(shù)

6)

svd(u,s,vt,CovArray)

//進(jìn)行奇異值分解求CovArray的特征值及特征向量，

//u為特征向量

7)

Y=matmul(ImageData,u)

1.3 基于HEVC的去空間冗余及編碼

1.3.1 HEVC框架

HEVC標(biāo)準(zhǔn)是ITU-T VCEG繼H.264之后所制定的新的視頻編碼標(biāo)準(zhǔn)，圍繞著現(xiàn)有的視頻編碼標(biāo)準(zhǔn)H.264進(jìn)行了優(yōu)化改良，編碼性能基本達(dá)到H.264的2倍。圖4描述了HEVC標(biāo)準(zhǔn)的主要編碼框架。

圖4 HEVC 編碼框架

1.3.2 HEVC中的量化系數(shù)(16 bit的擴(kuò)展)

HEVC使用YCbCr顏色空間，Y指亮度分量，Cb和Cr指色度分量。HEVC支持的采樣格式有4∶4∶4、4∶2∶2、4∶2∶0和4∶0∶0，適用于高光譜圖像的格式為4∶0∶0。HEVC定義的像素精度包括；8 bit、10 bit、12 bit，且已有開(kāi)源HEVC實(shí)現(xiàn)中最成熟的x265只實(shí)現(xiàn)了8 bit、10 bit兩種。高光譜圖像每個(gè)譜段僅有一個(gè)表示光譜反射率的分量，且常見(jiàn)的光譜儀拍攝到的數(shù)據(jù)多為16 bit。針對(duì)這種問(wèn)題，有兩種解決方法：1)將16 bit光譜數(shù)據(jù)進(jìn)行下采樣，轉(zhuǎn)化為新的8 bit光譜數(shù)據(jù)；2)擴(kuò)展x265在量化部分支持的位數(shù)。為保證壓縮質(zhì)量，本文選擇了第2種方案，對(duì)x265進(jìn)行了擴(kuò)展。

對(duì)DCT系數(shù)進(jìn)行量化的公式為：

Zij=(Wij×MF+2q-1)?q;q=21+[Qp/6]-lbN

(5)

其中：Zij表示量化后的殘差系數(shù)；Wij表示DCT后的殘差系數(shù)；Qp為量化參數(shù)；width為殘差矩陣的列數(shù)；MF為量化系數(shù)，其值可根據(jù)選定的Qp值查表得到。

在本文中，處理16 bit精度高光譜數(shù)據(jù)時(shí)所采用的量化系數(shù)MF矩陣經(jīng)計(jì)算后的具體值見(jiàn)式(6)。

(6)

HEVC編解碼器中，量化參數(shù)Qp和量化步長(zhǎng)Qstep的關(guān)系依如下規(guī)則確定：

1)量化步長(zhǎng)Qstep共有52個(gè)值(對(duì)于亮度編碼而言)。

2)量化參數(shù)Qp是量化步長(zhǎng)Qstep的序號(hào)，取值0～51。

3)Qp取最小值0 時(shí)，表示量化最精細(xì)；相反，Qp取最大值51時(shí)，表示量化是最粗糙的。

4)Qstep隨著Qp的增加而增加，Qp每增加6，Qstep增加1倍。

1.3.3 基于HEVC幀內(nèi)預(yù)測(cè)的去空間冗余及編碼

HEVC將原始圖像劃分為多個(gè)不同大小的編碼單元，根據(jù)圖像紋理復(fù)雜程度劃分為從64×64至8×8不等的尺寸。圖5為一幅高光譜圖像其中一個(gè)波段的子圖被HEVC幀內(nèi)預(yù)測(cè)算法編碼時(shí)的劃分情況。

HEVC利用預(yù)測(cè)+DCT消除空間冗余。根據(jù)紋理復(fù)雜度劃分不同大小的編碼單元，能夠提高DCT消除冗余的程度。

在對(duì)圖像進(jìn)行預(yù)測(cè)時(shí)，HEVC的預(yù)測(cè)方向多達(dá)34個(gè)，因此極大地提高了幀內(nèi)預(yù)測(cè)的準(zhǔn)確性，保證了高光譜圖像壓縮質(zhì)量。

從圖4可以看出，HEVC在完成預(yù)測(cè)、變換、量化后的操作為熵編碼，這也是生成二進(jìn)制碼流的關(guān)鍵步驟。HEVC采用基于內(nèi)容自適應(yīng)的二進(jìn)制算術(shù)編碼。

圖5 HEVC幀內(nèi)預(yù)測(cè)算法對(duì)一幅高光譜圖像的編碼單元?jiǎng)澐智闆r

2 原型系統(tǒng)實(shí)現(xiàn)與實(shí)驗(yàn)分析

采用Nvidia Jetson TX1 Developer Kit嵌入式硬件平臺(tái)對(duì)所設(shè)計(jì)算法進(jìn)行了實(shí)現(xiàn)。該硬件處理平臺(tái)核心硬件包含ARM架構(gòu)的A57型號(hào)4核處理器，Maxwell架構(gòu)含有256個(gè)cuda單元的Nvidia顯卡和4 GB LPDDR4內(nèi)存。

使用美國(guó)航天局提供的AVIRIS(Airborne Visible/InfRared Imaging Spectrometer)高光譜數(shù)據(jù)[18]作為測(cè)試進(jìn)行進(jìn)行性能驗(yàn)證。圖像場(chǎng)景為黃石公園，圖像大小為512×512，包含224個(gè)波段，每像素值占用16 bit，文件大小為112 MB。

作為對(duì)比算法的JPEG2000使用Kakadu 7.8[19]作為軟件實(shí)現(xiàn)。

對(duì)于高光譜圖像有損壓縮算法來(lái)說(shuō)，率失真是檢驗(yàn)其壓縮性能的重要標(biāo)準(zhǔn)之一。通常利用率失真來(lái)衡量壓縮后重建圖像與原始高光譜圖像的保真度，包括兩項(xiàng)指標(biāo)：均方差(Mean Squared Error，MSE)和峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)。

(7)

PSNR可定義為：

(8)

2.1 壓縮性能對(duì)比

為了達(dá)到和Kakadu 7.8相似的比特率，相應(yīng)地取x265的量化參數(shù)(Qp)為4、12、20、28。具體的x265軟件輸入?yún)?shù)見(jiàn)表1。

表1 x265參數(shù)信息

在不同碼率的情況下，本文所采用的算法與JPEG2000的性能比較如圖6所示?？梢钥闯?，本文所采用的算法在峰值信噪比意義下的性能表現(xiàn)均高于JPEG2000。

圖7(a)為在壓縮比為53∶1時(shí)的壓縮結(jié)果，文件大小為9.91 KB，PSNR=41.14 dB；圖7(b)為本文算法在壓縮比為53∶1時(shí)的壓縮結(jié)果，文件大小為9.8 KB，PSNR=43.52 dB。在肉眼直觀的判斷上，在圖7兩幅子圖左下角的地表細(xì)節(jié)部分，圖(b)明顯優(yōu)于圖(a)。

圖6 JPEG2000與本文算法的性能比較

圖7 JPEG2000與本文算法的壓縮效果對(duì)比

2.2 KTL并行計(jì)算時(shí)間分析

分別在CPU和GPU上對(duì)KLT進(jìn)行實(shí)驗(yàn)對(duì)比，結(jié)果如圖8所示，其中橫坐標(biāo)1～16分別表示512×512×10，512×512×20，…，512×512×160。從圖8可以看出，當(dāng)數(shù)據(jù)量較小，即高光譜圖像的譜段數(shù)少于70時(shí)，KLT在CPU的計(jì)算時(shí)間會(huì)少于GPU；當(dāng)高光譜圖像的譜段數(shù)超過(guò)70時(shí)，KLT在GPU的計(jì)算時(shí)間會(huì)少于CPU；隨著高光譜圖像的譜段數(shù)不斷增加，KLT在CPU上的計(jì)算時(shí)間呈線性增長(zhǎng)，而GPU的計(jì)算時(shí)間增加開(kāi)始變得緩慢，并慢慢趨于平穩(wěn)。實(shí)驗(yàn)結(jié)果表明，使用GPU進(jìn)行KLT明顯比CPU計(jì)算時(shí)間短，并且數(shù)據(jù)量越大，體現(xiàn)的效果越明顯，至多可縮短33%的時(shí)間。

圖8 KLT在CPU和GPU上計(jì)算時(shí)間對(duì)比

3 結(jié)語(yǔ)

本文使用GPU實(shí)現(xiàn)了KLT的并行化，并聯(lián)合使用KLT和HEVC提出了一種嵌入式實(shí)時(shí)高光譜壓縮算法。利用NVIDIA Jetson TX1平臺(tái)設(shè)計(jì)并實(shí)現(xiàn)了原型系統(tǒng)，基于真實(shí)數(shù)據(jù)集進(jìn)行了算法性能測(cè)試。實(shí)驗(yàn)結(jié)果表明，本文方案與DWT+JPEG2000方案相比，在低中高比特率下均具有巨大優(yōu)勢(shì)；同時(shí)通過(guò)縮短KLT的運(yùn)算時(shí)間，使系統(tǒng)滿足了實(shí)時(shí)要求。HEVC的關(guān)鍵幀壓縮算法具有較高的計(jì)算復(fù)雜度，對(duì)運(yùn)行平臺(tái)的性能有較高要求，使其能夠在低端設(shè)備上滿足實(shí)時(shí)性的優(yōu)化算法是一個(gè)未來(lái)需要研究的方向之一。