基于稀疏表示的高光譜數據壓縮算法

吳 倩 張 榮徐大衛

(中國科學技術大學電子工程與信息科學系 合肥 230027)

(中國科學院電磁空間信息重點實驗室 合肥 230027)

基于稀疏表示的高光譜數據壓縮算法

吳 倩 張 榮*徐大衛

(中國科學技術大學電子工程與信息科學系 合肥 230027)

(中國科學院電磁空間信息重點實驗室 合肥 230027)

如何降低高光譜圖像大規模數據的存儲和傳輸代價一直是學者們關心的問題。該文提出一種基于稀疏表示的高光譜數據壓縮算法，通過一種波段選擇算法構造訓練樣本集合，利用訓練得到的基函數字典對高光譜數據所有波段進行稀疏編碼，并對表示結果中非零元素的位置和數值進行量化和熵編碼，從而實現高光譜圖像壓縮。實驗結果表明該文算法與3維小波相比具有更好的非線性逼近性能，其率失真性能明顯優于3D-SPIHT，并且在光譜信息保留上具有巨大的優勢。

圖像處理；數據壓縮；高光譜遙感；稀疏表示

1 引言

高光譜遙感數據憑借其豐富的光譜信息而被廣泛應用于環境監測、地質調查、資源勘探等領域。然而，考慮到此類多維大規模數據高昂的存儲和傳輸代價，急需高效的壓縮技術對數據進行編碼。

國內外學者已經提出很多算法來解決這一問題。一般而言，高光譜遙感圖像壓縮算法可以分為兩類：基于預測的無損壓縮算法和基于變換的有損壓縮算法。基于預測的無損壓縮算法利用高光譜數據譜間相關性和空間相關性，通過特定的預測準則對當前像素進行估計并獲取預測誤差，如基于查找表的壓縮算法及其相關改進算法[1- 3]。最近的一些無損壓縮算法參考文獻[4,5]。在有損壓縮算法中，各種變換，如離散小波變換(DWT), KL變換等，被用來去除高光譜數據的譜間和空間冗余，然后對變換系數進行量化和熵編碼[6]。由于小波變換良好的數學特性及相關的基于濾波器和提升框架的快速算法的提出，基于小波變換的編碼技術在過去若干年處于統治地位[7- 10]。另外，由于PCA是最小均方誤差意義下的最優去相關變換，JPEG2000具有最好的碼率失真特性，PCA-JPEG2000[11]是目前最優的高光譜圖像壓縮方案，在低、中、高比特率下，都可以達到較為理想的壓縮效果。近幾年來，針對高光譜壓縮領域的研究主要集中于低復雜度算法及高效編碼方案的研究，同時有損壓縮算法對數據信息的提取能力也更受到重視。文獻[12]提出利用非線性PCA實現高光譜數據壓縮和目標檢測；文獻[13]將視頻編碼標準H.264/AVC應用于高光譜圖像，證明基于H.264/AVC方法設計低功率、實時的高光譜數據編碼器的可能性；更多研究參見文獻[14-17]。

然而，以上方法中所采用的固定的基函數對信號/圖像的稀疏化表示能力有限，制約壓縮比的進一步提高。與傳統的固定正交基相比，一個通過訓練得到的過完備的基函數字典可以更為稀疏地對信號進行表示。隨著稀疏表示理論的發展，一些基于稀疏表示的算法被提出用于壓縮自然圖像和人臉圖像[18]。考慮到高光譜數據極強的譜間相關性與結構相似性，將稀疏表示理論應用于高光譜圖像壓縮具有合理性。

本文提出了一種基于稀疏表示的高光譜數據壓縮算法，利用一個通過訓練得到的基函數字典對高光譜數據所有波段進行稀疏編碼，并對表示結果中非零元素的位置和數值進行量化和熵編碼，從而實現高光譜圖像壓縮。

2 稀疏表示理論

從某種角度來說，稀疏表示理論研究的是一個關于信號表示的數學模型。具體而言，對于信號x∈?n，其稀疏表示形式為

即信號x是字典中一系列原子{di}的線性組合，為各個原子對應的權重系數。

稀疏表示問題可以分為稀疏編碼和字典學習兩個部分。給定信號x和字典D，稀疏編碼就是尋找稀疏解ω的過程，它可以表示為如式(3)所示的優化問題：

式(3)所示l0范數最小化問題為NP難題，通常采用貪心算法近似求解該問題，如正交匹配追蹤(Orthogonal Matching Pursuit, OMP)[19]。與稀疏編碼不同，字典學習用于估計基函數字典D。給定一個包含M個信號的訓練樣本集合字典可以通過求解下面的優化問題獲得。

或

其中ε表示允許誤差。簡單地說，字典學習問題就是通過在字典D的基礎上，尋找每一個樣本xi最稀疏的表示系數ωi，以此獲得對訓練樣本信號的合理表示結果，并得到所需的字典。字典學習算法多種多樣，如奇異值分解法(K-SVD)[20]，遞歸最小二乘法(RLS)[21]，等等。

3 基于稀疏表示的高光譜圖像壓縮算法

高光譜遙感數據由同一地物對不同波長電磁波(可見光和近紅外波段)反射成像而得，各波段之間具有很強的結構相似性，如圖1所示。

根據稀疏表示相關理論，與所需要表示的信號相比，訓練樣本極具代表性，因而可以利用訓練所得到的基函數字典對信號進行高效、稀疏表示。另一方面，對高光譜遙感數據而言，不同波段圖像間存在很強的相關性與結構相似性，這使得利用高光譜數據某一波段數據作為樣本訓練得到基函數字典，然后對其他波段進行稀疏表示，并對稀疏表示系數進行熵編碼以實現數據壓縮。基于以上想法，本文提出了如下高光譜遙感數據壓縮方案，包括：字典學習、稀疏編碼和熵編碼3個部分，如圖2所示。

3.1 高光譜數據的稀疏表示方法

圖1 AVIRIS Cuprite場景不同波段圖像

圖2 基于稀疏表示的高光譜遙感數據壓縮方案流程

根據第2節字典學習理論，給定包含個L訓練樣本集合高光譜數據的字典可以通過式(6)所示的優化問題得到

圖3給出一個利用K-SVD[20]算法對AVIRIS (Airborne Visible/InfraRed Imaging Spectrometer)數據Cuprite場景進行學習得到的字典，圖中字典大小為D∈?256× 256。

圖3 字典學習結果

3.2 訓練波段選擇

本文將稀疏表示理論應用于高光譜數據壓縮，整個壓縮流程中，需要選擇一個合適的訓練波段作為訓練樣本來進行字典學習。一種簡單的做法就是隨機選擇一個波段作為訓練集。然而，不同的訓練波段對字典的表示性能會有不同的影響。對于某些波段，特別是大氣吸收作用明顯的波段，它們形似噪聲，與其他波段相關性極低，如圖4所示，如果選擇這些波段作為訓練波段，得到的字典的表示性能將不盡人意。

圖4 Cuprite場景“噪聲”波段

根據稀疏表示理論，當選擇某波段作為訓練樣本時，字典學習問題是通過最小化重建波段和原始波段之間的逼近誤差得到基函數字典。另一方面，稀疏編碼是在給定稀疏度或允許誤差的條件下根據訓練所得的字典找到稀疏解。不難理解，訓練波段與其需要用字典進行稀疏表示的波段之間相關性越高，訓練得到的字典將具有更好的表示性能。因此，可以通過選擇與其他波段相關程度最高的波段作為訓練波段。

相關系數是衡量兩個變量之間相關程度的指標。定義第k個波段和第k+s波段的相關系數為

其中σk, k+s為第k個波段和第k+s波段的協方差，σk和σk+s對應兩個波段的方差，如式(9)和式(10)所示，為第k個波段的均值。

如3.1節所述，高光譜數據的表示是通過對波段圖像分塊進行的，考慮到不同圖像塊之間的關系，利用塊相關系數對波段間的相關性重新定義：

表1 訓練波段選擇算法

值得注意的是，訓練波段選擇的目的是尋找一個可以接受的訓練集，以保證訓練樣本具有一定程度的代表性，從而確保對高光譜數據稀疏表示的高效性。事實上，所得到的訓練波段只是一個相對較優的選擇，但不一定是最優波段。

采用AVIRIS數據Cuprite, Jasper Ridge和Lunar Lake, 3個場景第11～42波段數據圖5所示的訓練波段算法的性能進行驗證。表2為幾種不同準則組合下波段選擇結果。其中C1和C2分別表示利用式(8)和式(11)計算相關系數。

由于高光譜數據相鄰波段數據相關程度極高，幾乎沒有區別，結合表2數據，除了通過波段選擇算法得到的訓練波段之外，我們選取了其他若干波段作為訓練波段，3個場景的率失真性能對比如圖5所示。

表2 波段選擇結果

圖5所示結果表明，經過不同訓練波段學習得到的字典的率失真性能不同，其中C2+Mean的方式得到結果更為出色，并且計算簡單。因此，我們建議采用C2+最大化相關系數均值的方式選擇訓練波段。

3.3 壓縮方案

圖5 不同訓練波段情況下率失真性能比較

采用與前文相同的定義，假設稀疏編碼后得到的稀疏系數矩陣為由于ω具有很高的稀疏性，即只有很少的元素為非零元素。那么，對于每一個數據塊的稀疏解可用如下向量表示，其中分別表示第i個非零元素的位置索引值和數值，并且，由于對和的量化熵編碼有兩種方式：一種是按照中元素大小位置將和重排，對重排后的進行差分編碼，對量化熵編碼；另一種則是按照非零元素的賦值大小，即按中元素大小順序對將和重排，對重排后的進行差分編碼，對進行量化熵編碼。對全部稀疏編碼波段而言，最終的稀疏解為表示每個波段的稀疏解。

這里采用第1種方式對高光譜數據稀疏表示結果進行量化編碼。對非零元素字典序數由于D∈?n× m, pi可能取到的最大值為m，那么，不考慮差分的影響，編碼一個數據塊所需的非零元素位置序號所需比特數為假設記錄一個非零元素值所需的比特數為1，那么，整個壓縮的比特率R為

4 實驗與分析

以下實驗采用1997年AVIRIS的Cuprite，Jasper Ridge和Lunar Lake 3個場景的數據。AVIRIS圖像包含224個波段，其寬度是614。為了方便起見，每個波段圖像從左上角開始被切成512×512的子圖像。本文實驗中采用K-SVD訓練字典，以下實驗中所用字典D∈?256× 1024，即字典包含1024個原子，波段圖像分塊大小為16×16。

4.1 壓縮性能的評價指標

重構的圖像的信噪比(SNR)通常被用作評價單一波段圖像的壓縮性能的標準之一。假設高光譜圖像重建波段為第k個波段的信噪比定義為

對于高光譜數據，我們采用各重建波段的平均信噪比(ASNR)衡量整個算法的壓縮性能：其中t為訓練波段。

高光譜遙感技術最重要的應用之一是分類，其豐富的光譜信息可以用于識別不同類型的地物。從這個角度來看，對光譜信息的高保真重建是高光譜數據壓縮需要滿足的要求之一。光譜角是進行光譜匹配或高光譜數據分類的一個準則，它可以反映兩條光譜曲線的相似程度。對兩條光譜曲線其光譜角SA (Spectral Angle)定義為兩個向量的夾角，如式(15)所示。

其中，＜·＞表示兩個向量的內積。

這里用平均光譜角ASA(Average Spectral Angle)作為評價標準，定義如式(16)，以此來衡量壓縮算法對高光譜數據光譜信息的保真程度。

其中SAi, j表示空間位置為(i, j)的重建像元與原始像元之間的光譜角。

4.2 實驗結果和討論

4.2.1逼近性能比較 利用訓練所得字典對高光譜數據進行稀疏表示，并與3維小波(3D-DWT)的逼近性能進行對比，結果如圖6所示。圖6中曲線為不同有效系數個數百分比(即保留的變換域系數的個數與全部系數個數之比)情況下對應重建圖像的SNR曲線。其中3D-DWT采用9/7雙正交小波基，分別進行3級分解(LEVEL 3)與5級分解(LEVEL 5)。

圖6結果表明，在保留25%以內的有效系數時，稀疏表示方法的逼近性能優于3維小波的逼近性能，特別地，在僅保留少量系數的情況下，稀疏表示方法的逼近性能更為突出，這從側面反映了字典學習表示算法在低碼率壓縮方面的優勢。

圖6 非線性逼近性能對比

4.2.2率失真性能比較 本節對不同壓縮碼率下稀疏表示壓縮方案與3D-SPIHT壓縮算法的壓縮性能進行比較。圖7顯示了具體的壓縮碼率對比結果，圖8為0.4 bpp情況下不同算法重建圖像視覺效果對比。其中3D-SPIHT采用開源軟件QccPack實現，采用9/7小波，對應的兩條曲線中LEVEL 3和LEVEL 5分別表示算法中的3D-SPIHT所使用的小波分別進行3級分解和5級分解。本文算法里采用均勻量化的方式量化表示系數，并對其進行算術編碼。實驗結果表明本文算法明顯優于3D-SPIHT，且對圖像的結構信息具有更好的保真性能。

4.2.3光譜信息保留性能 為了評價壓縮算法對高光譜數據光譜信息的影響大小，下面對不同比特率下重建數據與原始數據的光譜角進行對比，結果如表3所示。表3結果表明，基于稀疏表示的高光譜壓縮算法在對原始光譜信息的保留方面具有巨大優勢。圖9給出了重建數據與原始數據光譜曲線對比。很明顯，本文算法重建數據更貼近原始數據，對高光譜圖像光譜信息的保留能力更強。

5 結束語

本文利用高光譜數據各波段的相似性，結合稀疏表示理論，提出了一種基于稀疏表示的高光譜數據壓縮方案。實驗結果表明，本文方案與以經典的3D-SPIHT算法相比，在中低比特率下具有巨大優勢，并且本文算法可以更好地保留高光譜數據的光譜信息，對壓縮后處理工作具有一定意義。由于稀疏表示算法具有較高的計算復雜度，其快速算法，特別是適用于高光譜數據的優化算法，是一個值得研究的方向。同時，與自然圖像相比，高光譜圖像具有更為復雜的結構特征，基于多尺度字典的壓縮算法是未來需要研究的工作之一。

圖7 率失真性能比較

圖 8 不同算法重建圖像視覺效果對比，各子圖為0.4 bpp情況下不同算法Band 120的重建圖像

表3 光譜角對比(°)

圖9 重建圖像光譜曲線對比

[1] Huang B and Sriraja Y. Lossless compression of hyperspectral imagery via lookup tables with predictor selection[J]. Remote Sensing, 2006, 63650L-63650L-8.

[2] Mielikainen J. Lossless compression of hyperspectral images using lookup tables[J]. IEEE Signal Processing Letters, 2006, 13(3): 157-160.

[3] Mielikainen J and Toivanen P. Lossless compression of hyperspectral images using a quantized index to lookup tables[J]. IEEE Geoscience and Remote Sensing Letters, 2008, 5(3): 474-478.

[4] Mielikainen J and Huang B. Lossless compression of hyperspectral images using clustered linear prediction with adaptive prediction length[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(6): 1118-1121.

[5] Cheng K and Dill J. Hyperspectral images lossless compression using the 3D binary EZW algorithm[C]. SPIE Electronic Imaging, International Society for Optics and Photonics, 2013: 865515-865515-8.

[6] Penna B, Tillo T, Magli E, et al.. Transform coding techniques for lossy hyperspectral data compression[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(5): 1408-1421.

[7] Tang X and Pearlman W A. Three-dimensional Wavelet-Based Compression of Hyperspectral Images[M]. Hyperspectral Data Compression, Springer US, 2006: 273-308.

[8] Rucker J T, Fowler J E, and Younan N H. JPEG2000 coding strategies for hyperspectral data[C]. Proceedings of IEEE International Geoscience and Remote Sensing Symposium, Seoul, Korea, 2005: 1-4.

[9] Penna B, Tillo T, Magli E, et al.. Progressive 3-D coding of hyperspectral images based on JPEG 2000[J]. IEEE Geoscience and Remote Sensing Letters, 2006, 3(1): 125-129. [10] Chang Chein-I. Hyperspectral Data Exploitation: Theory and Applications[M]. New York: John Wiley, 2007: 379-407. [11] Du Q and Fowler J E. Hyperspectral image compression using JPEG2000 and principal component analysis[J]. IEEE Geoscience and Remote Sensing Letters, 2007, 4(2): 201-205. [12] Du Q, Wei W, Ma B, et al.. Hyperspectral image compression and target detection using nonlinear principal component analysis[J]. SPIE International Society for Optics and Photonics, 2013: 88710S-88710S-6.

[13] Santos L, López S, Callico G M, et al.. Performance evaluation of the H. 264/AVC video coding standard for lossy hyperspectral image compression[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012, 5(2): 451-461.

[14] Wang C, Miao Z, Feng W, et al.. An optimized hybrid encode based compression algorithm for hyperspectral image[C]. International Conference on Optical Instruments and Technology (OIT2013), Beijing, China, 2013: 90451V-90451V-7.

[15] Cheng K and Dill J. Hyperspectral images lossless compression using the 3D binary EZW algorithm[J]. International Society for Optics and Photonics/SPIE Electronic Imaging, 2013: 865515-865515-8.

[16] Pan X, Liu R, and Lü X. Low-complexity compression method for hyperspectral images based on distributed source coding[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(2): 224-227.

[17] Noor N R M and Vladimirova T. Integer KLT Design Space Exploration for Hyperspectral Satellite Image Compression [M]. Convergence and Hybrid Information Technology, Berlin Heidelberg Springer, 2011: 661-668.

[18] Bryt O and Elad M. Compression of facial images using the K-SVD algorithm[J]. Journal of Visual Communication and Image Representation, 2008, 19(4): 270-282.

[19] Tropp J A and Gilbert A C. Signal recovery from random measurements via orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655-4666.

[20] Aharon M, Elad M, and Bruckstein A. K-svd: an algorithm for designing overcomplete dictionaries for sparse representation[J]. IEEE Transactions on Signal Processing, 2006, 54(11): 4311-4322.

[21] Skretting K and Engan K. Recursive least squares dictionary learning algorithm[J]. IEEE Transactions on Signal Processing, 2010, 58(4): 2121-2130.

吳 倩： 女，1988年生，博士生，研究方向為遙感圖像處理、數據壓縮.

張 榮： 女，1968年生，副教授，研究方向為數字圖像處理、遙感圖像處理、數據壓縮.

徐大衛： 男，1990年生，碩士生，研究方向為遙感圖像處理.

Hyperspectral Data Compression Based on Sparse Representation

Wu Qian Zhang Rong Xu Da-wei
(Department of Electronic Engineering and Information Science, University of Science and Technology of China, Hefei 230027, China)
(Key Laboratory of Electromagnetic Space Information, Chinese Academy of Science, Hefei 230027, China)

How to reduce the storage and transmission cost of mass hyperspectral data is concerned with growing interest. This paper proposes a hyperspectral data compression algorithm using sparse representation. First, a training sample set is constructed with a band selection algorithm, and then all hyperspectral bands are coded sparsely using a basis function dictionary learned from the training set. Finally, the position indices and values of the non-zero elements are entropy coded to finish the compression. Experimental results reveal that the proposal algorithm achieves better nonlinear approximation performance than 3D-DWT and outperforms 3D-SPIHT. Besides, the algorithm has better performance in spectral information preservation.

Image processing; Data compression; Hyperspectral remote sensing; Sparse representation

TP751.1

A

1009-5896(2015)01-0078-07

10.11999/JEIT140214

2014-02-19收到，2014-05-20改回

國家自然科學基金(61172154)，國家973計劃項目(2010CB731904)，國家自然科學基金重點項目(61331020)和中國科學院光電研究院雛鷹計劃資助課題

*通信作者：張榮 zrong@ustc.edu.cn