利用格型遞歸最小二乘濾波器組的高光譜圖像壓縮

鄭 鐵，薛長斌，宋金偉

（1.中國科學(xué)院 國家空間科學(xué)中心，北京100190；2.中國科學(xué)院大學(xué)，北京100190；3.國家電網(wǎng)有限公司大數(shù)據(jù)中心，北京100031）

1 引言

高光譜圖像被廣泛應(yīng)用于地質(zhì)勘探、環(huán)境分析、農(nóng)業(yè)監(jiān)測等領(lǐng)域。隨高光譜圖像分辨率的提升，圖像數(shù)據(jù)量急劇增加，給傳輸帶寬和存儲容量帶來巨大挑戰(zhàn)。因此，高效的高光譜圖像壓縮技術(shù)必不可少［1-5］。高光譜圖像的壓縮方法分為有損壓縮和無損壓縮。有損壓縮雖然能達到很高的壓縮比，但會破壞數(shù)據(jù)的細節(jié)信息，這對需要精確分析的計算場景是難以接受的，只有無損壓縮能保證圖像信息的完整性［6-7］。

高光譜圖像同時具備空間相關(guān)性和譜間相關(guān)性。傳統(tǒng)的二維圖像壓縮技術(shù)只針對空間相關(guān)性，無法利用譜間相關(guān)性提升壓縮比［8-10］。目前，常見的高光譜圖像壓縮方法分為變換類、矢量量化類和預(yù)測類。其中變換類算法是先將圖像從空域映射到變換域，使用較少碼字表示圖像的主要信息，隨后再對細節(jié)量化編碼，因此常用于有損壓縮［11］。矢量量化類算法的計算復(fù)雜度會隨維數(shù)增加呈現(xiàn)指數(shù)級的增長，龐大的計算量限制了該類算法的實際應(yīng)用［12］。預(yù)測類算法先使用預(yù)測器計算待測像素點估計值與實際值的差異，去除圖像的相關(guān)性，再針對預(yù)測誤差做壓縮編碼。這類算法在獲取良好壓縮結(jié)果的同時保持較低的算法復(fù)雜度，逐漸成為國內(nèi)外研究熱點［13］。

Mielikainen等人提出基于查找表（Lookup Table，LUT）的預(yù)測方案［14］，該方案將前一波段中與當(dāng)前像素點具有同位坐標的像素作為索引，使用查找表返回預(yù)測結(jié)果。H.Bormin等人改進LUT的預(yù)測方案，提出局部平均譜間尺度變換的方法（Locally Averaged Inter-Band Scaling-LUT，LAIS-LUT）［15］。該方法先計算待測像素點的LAIS估計值，隨后在兩個查找表中索引與LAIS估計值最接近的值并作為預(yù)測結(jié)果?；贚UT的壓縮方法雖然具有良好的預(yù)測效果，然而它消耗大量的內(nèi)存資源，不適用于星載圖像壓縮系統(tǒng)。

聚類差分脈沖編碼調(diào)制（Clustered Differ?ence Pulse Code Modulation，C-DPCM）是一種將聚類和線性濾波器相結(jié)合的預(yù)測方案［16］。該方法首先利用k-means聚類算法將光譜矢量做聚類處理，提升同類矢量間的相似度。隨后，采用預(yù)測長度固定的線性濾波器做預(yù)測壓縮。具有自適應(yīng)預(yù)測長度的聚類算法（C-DPCM with Adaptive Prediction Length，C-DPCM-APL）是C-DPCM算法的擴展，它通過搜尋線性濾波器輸入的最優(yōu)預(yù)測長度獲得良好的壓縮效果［17］。但該方案的計算復(fù)雜度較高、消耗資源較大，更適合被應(yīng)用于離線壓縮場景。

空間數(shù)據(jù)系統(tǒng)協(xié)商委員會（Consultative Committee for Space Data Systems，CCSDS）提出基于最小均方誤差濾波器的快速無損（Fast lossless，F(xiàn)L）算法［18］，該算法利用先前的波段預(yù)測當(dāng)前像素，并在線性預(yù)測過程中自適應(yīng)改變預(yù)測系數(shù)。三階譜間預(yù)測器（The third-order interband predictor，IP3）［19］由兩階段預(yù)測器和熵編碼組成，通過譜間預(yù)測器和反向像素搜索（BPS）方案計算最終預(yù)測值。宋金偉等人提出基于遞歸最小二乘法（Recursive Least Squares，RLS）的壓縮方案，該方案先將待測像素點做空間濾波，再利用固定長度的遞歸最小二乘做譜間預(yù)測［20］。由于遞歸最小二乘濾波器較最小均方誤差濾波器的收斂速度快，因此其壓縮結(jié)果要優(yōu)于FL算法。高放和郭樹旭提出自適應(yīng)預(yù)測波段的傳統(tǒng)最小二乘壓縮算法（Cnventional Recursive Least-Squares Predictor with Adaptive prediction bands，ACRLS）［21］。該文獻指出不同待測波段的最優(yōu)預(yù)測長度是不盡相同的，并通過窮舉法搜索每個待測點RLS濾波器的最優(yōu)預(yù)測長度，獲得良好的壓縮效果。但該算法的計算復(fù)雜度過高，消耗大量的計算時間。宋金偉等人提出的快速RLS（Fast Adaptive-Length-Prediction RLS，F(xiàn)ast-RLS-ALP）預(yù)測算法［22］，相較于傳統(tǒng)RLS算法，該算法充分利用已知波段的濾波器做前向預(yù)測，降低計算復(fù)雜度。但該算法的空間預(yù)測方案比較簡單，且譜帶間濾波器只使用固定長度的預(yù)測方案，因此壓縮結(jié)果并不突出。

本文提出一種基于格型遞歸最小二乘濾波器組的高光譜圖像無損壓縮方案（RLS Lattice Filter Group，LFG-RLS）。該壓縮方案采用的格型濾波器組充分考慮最優(yōu)譜帶預(yù)測長度對預(yù)測效果的影響，合理地為每個譜帶選擇最優(yōu)濾波器，并利用其鏈式序列更新的特點，對濾波器組分類、簡化，大幅度降低計算時間。

2 壓縮方案介紹

本文提出的壓縮方案分為預(yù)測器和編碼器，其中預(yù)測器包含譜帶內(nèi)預(yù)測和譜間預(yù)測兩個階段。對于高光譜圖像的每個譜帶，首先使用單邊高斯預(yù)測器做譜帶內(nèi)預(yù)測，得到譜帶內(nèi)預(yù)測誤差。其次，采用格型濾波器組篩選每個譜帶的最優(yōu)濾波器，使壓縮結(jié)果接近最優(yōu)，并根據(jù)格型濾波器組鏈式更新的特點，對最優(yōu)濾波器做進一步篩選，大幅度降低計算時間。最后，利用算術(shù)編碼器對預(yù)測誤差做壓縮編碼，得到壓縮碼流。本文的壓縮方案結(jié)構(gòu)圖如圖1所示。

圖1 壓縮方案結(jié)構(gòu)圖Fig.1 Algorithm structure chart

2.1 譜帶內(nèi)預(yù)測

高光譜圖像可被視為一組三維數(shù)據(jù)，設(shè)其尺寸為W×H×N（行×列×波段），本文使用sz(x，y)表示第z波段中第x行第y列像素點的灰度值，或使用sz(n)表示z波段中第n個像素點的灰度值，其中n=W×y+x。并將具有相同y坐標的平面稱為線譜平面（Line-Spectral Plate，LSP），如圖2所示。

圖2 線譜平面Fig.2 Line-Spectral Plate

在譜帶內(nèi)預(yù)測階段，通過對待測像素點的上下文像素做預(yù)測運算，計算譜帶內(nèi)預(yù)測誤差，從而去除高光譜圖像的空間相關(guān)性，其中待測像素點的上下文窗口如圖3所示。ACRLS［21］和Fast-RLS-ALP［22］算法分別采用24鄰域和4鄰域均值的預(yù)測方式去除空間冗余信息。

圖3 待測像素點的上下文窗口Fig.3 Context window of current pixel

因待測像素點與上下文像素之間的相關(guān)性隨距離的增加而減小，高斯預(yù)測器被用于譜帶內(nèi)預(yù)測。為保證預(yù)測過程的因果性，即使用已知像素點預(yù)測當(dāng)前待測點，采用的單邊高斯預(yù)測器如式（1）所示：

其中：

則 當(dāng) 前 像 素 點 的 估 計 值s?z(x，y)如 式（3）所示：

z波段第n個像素點的譜帶內(nèi)預(yù)測誤差d z(n)：

對于首個波段的預(yù)測誤差d z(n)直接使用編碼器壓縮編碼，其余波段則使用格型濾波器組去除譜間冗余信息后再編碼。

2.2 譜帶間預(yù)測

高光譜圖像在光譜維度上具有較強的相關(guān)性，但受噪聲影響，各相鄰波段譜間相關(guān)性的強弱不同。因此，RLS濾波器的預(yù)測效果與預(yù)測長度（待測像素點濾波器輸入向量的長度）并不總是呈現(xiàn)正相關(guān)的關(guān)系，對于部分譜帶，較短的預(yù)測長度可以獲得良好的預(yù)測效果［17，21］，合理地為每個譜帶的濾波器選擇預(yù)測長度，是提升壓縮算法性能的關(guān)鍵。

2.2.1 格型濾波器組

對于波段按行交叉格式的高光譜圖像，在譜帶間預(yù)測過程中，z波段第n個像素點共使用z-1個不同預(yù)測長度的濾波器F z，i(n)(i={1，2，…，z-1})，其中下標“i”被用于標識不同預(yù)測長度的濾波器。則F z，i(n)的輸入向量xz，i(n)=[d z-i(n)，…，d z-1(n)]T，對應(yīng)的權(quán)重向量w z，i(n)=[wz-i(n)，…，w z-1(n)]T，期 望 信 號為d z(n)。則F z，i(n)的計算過程如公式（5）～公式（13）所示。

首先計算Fz，i(n)的譜帶間預(yù)測殘差：

預(yù)測誤差方差的累計值vz，i(n)：

設(shè)：

其 中：cz，i(n)，uz，i(n)分 別 為 濾 波 器 的 投 影 向 量、投影值。

濾波器的增益為：

則權(quán)重向量w z，i(n)的遞歸更新如式（10）所示：

為確保濾波器的收斂性及屏蔽一些異常的像素值，當(dāng)濾波器階數(shù)≤3時，投影向量cz，i(n)及投影值uz，i(n)采用傳統(tǒng)RLS濾波方式更新，如式（11）～式（13）所示。

其中：根據(jù)文獻［21］RLS濾波器敏感性測試結(jié)果初 始 化 遺 忘 因 子 為λ=0.9995；w z，i(0)=[0，0，…，0]T。依 據(jù) 上 述 更 新 過 程 可 知：F z-i+1，1(n)，…，F(xiàn) z，i(n)是依次遞歸更新的，這組具有前向迭代更新關(guān)系的濾波器稱為鏈式序列。因高光譜圖像共有N個波段，首個波段的譜帶內(nèi)預(yù)測誤差直接使用編碼器壓縮編碼，本文提出的格型濾波器組中共包含N-1條鏈式序列，如圖4所示，為簡潔表示，省略每個濾波器的部分符號(n)。圖中每一列表示同一待測點不同預(yù)測長度的濾波器；箭頭表示同一條鏈式序列內(nèi)濾波器的迭代更新方向，將用于迭代更新F z，i(n)的鏈式序列記作集合Ωz，i(n)：

圖4 格型濾波器組的鏈式序列Fig.4 Chain sequence for filters

格型濾波器組在預(yù)測過程中，v z，i(n)被用于衡量同一待測點不同預(yù)測長度濾波器的預(yù)測精度。預(yù)測誤差被視為服從均值為0的高斯分布，其信息熵H與方差σ2呈現(xiàn)正相關(guān)關(guān)系，如式（15）所示：

因此，在波段z第n個像素點的不同預(yù)測長度濾波器中，具有最小預(yù)測誤差方差累計值v z，min(n)的濾波器稱為最優(yōu)濾波器F z，opt(n)，其中v z，min(n)=min{v z，1，…，vz，z-1(n)}。最優(yōu)濾 波 器的預(yù)測長度被稱為最優(yōu)預(yù)測長度，記作L z(n)。

若直接使用格型濾波器組順次計算每個待測點不同預(yù)測長度的濾波器，雖然能獲得最佳的預(yù)測效果，但耗費過長的計算時間。在此，通過對濾波器收斂性的考察，優(yōu)化了格型濾波器組的更新方式，為衡量z波段第n個與第n-1個待測點權(quán)值的變化情況，定義相鄰濾波器權(quán)值的距離（Distance of Weights，DOW）

以AVIRIS 2006高光譜圖像集中Yellow?Stone 0校準圖像為例，濾波器F223，222(n)的DOW隨迭代次數(shù)的變化趨勢如圖5所示，隨迭代次數(shù)的增加，各個濾波器具有明顯的收斂趨勢。因此，同一譜帶不同預(yù)測長度濾波器參數(shù)v z，i(n)的排布次序趨于穩(wěn)定，最優(yōu)濾波器搜索階段可以被提前終止，后續(xù)譜帶的最優(yōu)濾波器被認為是相同的。

圖5 不同LSP中相鄰濾波器權(quán)值的距離Fig.5 DOWsof filtersfor different LSPin calibrated image

引入限定迭代次數(shù)的閾值T，當(dāng)n=T時，先計算波段z第T個待測點的最優(yōu)濾波器F z，opt(T)，隨 后 更 新 其 鏈 式 序 列 集 合Ωz，opt(T)。為統(tǒng)計當(dāng)前譜帶所有待測點最優(yōu)濾波器的鏈式序列，使用最優(yōu)濾波器集合Ω表示每個待測點Ωz，opt(T)(1＜z≤N)的并集。當(dāng)n＞T時，最優(yōu)譜帶搜索階段將被終止，進入快速預(yù)測階段，只有位于Ω中的濾波器需要被持續(xù)更新。為直觀地描述格型濾波器組的預(yù)測過程，算法流程圖如圖6所示。

圖6 格型濾波器組的預(yù)測過程Fig.6 Prediction process for Lattice Filters

在上述格型濾波器組的更新過程中，當(dāng)n＞T時，位于Ω中的濾波器需要被持續(xù)更新。因此，越多的濾波器位于同一條鏈式序列中，待更新濾波器的數(shù)量越少，消耗的計算時間越短。

2.2.2 提升策略

以校準后Yellow Stone 0圖像的400th和500thLSP為例，各譜帶的最優(yōu)預(yù)測長度如圖7所示，除個別受噪聲影響的波段外，大多數(shù)譜帶的最優(yōu)預(yù)測長度隨z的增大而增大。因此，可以經(jīng)篩選分類，使更多的待更新濾波器處于同一鏈式序列。本文提出最長更新規(guī)則：若多個濾波器的最優(yōu)預(yù)測長度相近時，應(yīng)將具有最長鏈式序列的濾波器視為最優(yōu)濾波器。為簡化分類操作，若z波段第n個待測點的最優(yōu)濾波器預(yù)測長度L z(n)≥時，則將最長鏈式序列中的濾波器視為最優(yōu)濾波器；否則vz，i(n)最小的濾波器被用于迭代計算。

圖7 不同譜帶的最優(yōu)預(yù)測長度Fig.7 Optimal prediction length of different band

此外，在文獻［21］中，高放等人發(fā)現(xiàn)對于同一成像系統(tǒng)拍攝的高光譜圖像集，其最優(yōu)預(yù)測長度高度一致，并建議將首個圖像的最優(yōu)預(yù)測長度直接應(yīng)用于預(yù)測該場景的其余圖像。該加速方案同樣應(yīng)用于本文算法，記為RLS Speed Lattice Filter Group（SLFG-RLS）。

2.3 算術(shù)編碼

本文采用算術(shù)編碼器做壓縮編碼。因最終預(yù)測誤差數(shù)據(jù)是用16位表示的，故算術(shù)編碼的碼表需要使用65536個符號。根據(jù)預(yù)測殘差的概率分布，大多數(shù)符號沒有被使用，故采用自適應(yīng)碼表。初始的碼表包含0和ESC兩個符號。當(dāng)需要對一個新符號編碼時，編碼器先使用ESC的概率對其編碼，隨后再更新碼表便于后續(xù)編碼。

3 實驗結(jié)果與分析

本文壓縮方案中的所有程序均在實驗平臺（Intel i7-7700k CPU 4.20GHZ/16GBRAM）中以C語言進行編寫、測試。CCSDS的多光譜和高光譜數(shù)據(jù)壓縮工作組提供的AVIRIS 2006高光譜圖像數(shù)據(jù)集作為本文算法的測試數(shù)據(jù)源。

測試圖像集包括5個16位校準和5個16位未校準的場景，每個場景包含512行，224個譜帶的數(shù)據(jù)，其中校準與非校準場景每行分別有677和680個像素點。在此，將場景名稱的首字母及其序列號用作縮寫，例如：Yellowstone scene 0記作YS_SC0。

3.1 高斯預(yù)測器的參數(shù)選擇

在譜帶內(nèi)預(yù)測階段，為便于計算，本文采用整數(shù)型的單邊高斯預(yù)測器。本文將均值濾波器及不同參數(shù)的高斯預(yù)測器做測試與比較。以YS_SC0場景為例，實驗結(jié)果如表1所示。

表1 不同參數(shù)的高斯預(yù)測器Tab.1 Different parameters in Gaussian predictor

在窗口半徑相同時，對于每個譜帶高斯預(yù)測器比均值預(yù)測器多消耗十余毫秒的計算時間，但高斯預(yù)測器能獲得更佳的壓縮效果。對于高斯預(yù)測器，當(dāng)σ≥2時，能獲得較好的壓縮效果；當(dāng)窗口半徑r≥4時，預(yù)測器的壓縮效率沒有明顯提高，但時間復(fù)雜度卻驟增?？紤]到算法的計算復(fù)雜度和實時性，最終選用σ=2，r=4作為單邊高斯預(yù)測器的初始化參數(shù)。

3.2 迭代閾值的選擇

隨迭代次數(shù)的增大，濾波器的最優(yōu)預(yù)測長度趨于穩(wěn)定。以校準后YS_SC0場景中209，140，94，49波段的最優(yōu)濾波器為例，圖8展示不同LSP的最優(yōu)預(yù)測長度L z(n)。

圖8 不同LSP的最優(yōu)預(yù)測長度Fig.8 Optimal prediction length for different LSP

當(dāng)?shù)螖?shù)的閾值T≥50 LSPs時，最優(yōu)預(yù)測長度波動很小。根據(jù)公式（15）和公式（16），隨濾波器迭代訓(xùn)練次數(shù)的增加，收斂趨勢愈加明顯，同一譜帶不同預(yù)測長度濾波器參數(shù)vz，i(n)的排布次序及濾波器最優(yōu)預(yù)測長度波動漸小，趨于穩(wěn)定。為量化描述迭代閾值對算法的壓縮結(jié)果、計算時間的影響，以校準后YS_SC0場景為例，表2展現(xiàn)格型濾波器組在使用提升策略前的壓縮結(jié)果、計算時間與T的關(guān)系。

從表2中可知，當(dāng)T≥50 LSPs時，算法的壓縮結(jié)果趨于固定，但計算時間急劇上升。考慮到實際壓縮結(jié)果和時間復(fù)雜度的雙重影響，本文算法將參數(shù)T初始化為50 LSPs，從而保證良好的壓縮效果以及較低的時間復(fù)雜度。

表2 不同參數(shù)下的壓縮性能Tab.2 Compression performance under different pa?rameters

3.3 壓縮性能的評估與分析

為衡量本文算法的壓縮性能，將其與現(xiàn)有的 算 法（FL，IP3，RLS，ACRLS，and Fast-RLS-ALP）從壓縮結(jié)果與計算時間兩方面作對比。

表3展示各類算法的壓縮結(jié)果。其中，ACRLS算法因采用窮舉法查找每個波段的最佳預(yù)測長度，是已知同類型壓縮結(jié)果最佳的壓縮算法。本文提出的LFG-RLS算法采用格型濾波器組獲取圖像譜帶的最優(yōu)預(yù)測長度，并結(jié)合濾波器的收斂性，合理地選擇迭代閾值，獲取良好的壓縮效果。在校準、未校準場景中的平均壓縮結(jié)果分別是是3.34 bits/pixel和5.61 bits/pixel。該壓縮效果與ACRLS非常接近，并優(yōu)于其他壓縮方案。此外，SLFG-RLS的壓縮結(jié)果證實，當(dāng)首個場景的預(yù)測長度應(yīng)用于其他場景時，能獲得相似的壓縮效果。

表3 不同算法的壓縮結(jié)果Tab.3 Compression results of different algorithms （bits·pixel-1）

各算法的平均計算時間如表4所示，F(xiàn)L，IP3，RLS和Fast-RLS-ALP算法采用固定預(yù)測長度，預(yù)測方式較為簡單，計算時間較短。ACRLS算法采用窮盡式的搜索方案計算所有待測點的最優(yōu)預(yù)測長度，消耗較長的計算時間。

表4 各種算法的平均計算時間Tab.4 Average computation time of various algorithms （s）

本文提出的LFG-RLS算法，對于校準、未校準圖像分別僅需1050 s，1064 s，相較于使用窮舉法搜索最優(yōu)譜帶預(yù)測長度的ACRLS算法降低10余倍的運算時間。主要原因是本文合理地限定濾波器的迭代次數(shù)，并根據(jù)格型濾波器組鏈式更新的特點，對最優(yōu)濾波器做分類與簡化，大幅度減少待更新濾波器的數(shù)量，降低快速預(yù)測階段的計算時間。

此外，譜帶間預(yù)測器根據(jù)迭代閾值T可分為最優(yōu)譜帶搜索階段和快速預(yù)測階段，對于校準圖像，LFG-RLS壓縮方案兩階段平均耗時分別為767.8 s，262.9 s，最優(yōu)譜帶搜索階段壓縮的數(shù)據(jù)量約占總數(shù)據(jù)量的10%，但消耗總計算時間的73.1%。因此，若高光譜成像儀采集的圖像尺寸越長或圖像數(shù)量越多，最優(yōu)譜帶搜索階段消耗計算時間所占比例越小，壓縮方案所需總運算時間相對越短，能充分發(fā)揮算法的優(yōu)越性。在SLFGRLS壓縮方案中，系列校準、未校準圖像分別是同一系統(tǒng)采集的多幅圖像，將首個場景的預(yù)測結(jié)果推廣后，對校準、未校準圖像的平均計算時間分別降低至441 s和447 s。

為驗證本文壓縮算法的無損性，將壓縮后的數(shù)據(jù)流做解壓縮處理，得到恢復(fù)圖像。將原始圖像與恢復(fù)圖像相比較，結(jié)果顯示二者完全相同，表明高光譜圖像信息沒有任何損失。

4 結(jié)論

為充分挖掘高光譜圖像的空間相關(guān)性和譜間相關(guān)性，本文提出利用格型遞歸最小二乘法濾波器組的高光譜圖像無損壓縮方案。該方案中充分考慮最優(yōu)濾波器對預(yù)測效果的影響，采用格型濾波器組計算每個譜帶的最優(yōu)濾波器，并且根據(jù)其鏈式更新的特點，對篩選過程做進一步的分類與簡化，大幅度降低壓縮方案的計算復(fù)雜度。實驗結(jié)果表明，以CCSDS提供的AVIRIS高光譜圖像為測試數(shù)據(jù)集，本文的壓縮方案具有良好的壓縮性能，對16位校準圖像、未校準圖像分別取得3.34 bits/pixel和5.61 bits/pixel的壓縮結(jié)果，與同類算法的最佳壓縮結(jié)果相似，且其平均計算時間遠低于其余具有相似壓縮結(jié)果的算法，具有較強的競爭力，為高光譜圖像無損壓縮提供了一種有效的解決方案。未來將針對壓縮方案的并行化處理做研究與探討，并開發(fā)可行的硬件實現(xiàn)方案。