基于張量分解的光譜圖像壓縮感知重構 *

趙梓淵 ，唐意東 ，黃樹彩

（1. 空軍工程大學 信息與導航學院，陜西 西安 710077；2. 中國人民解放軍95607 部隊，四川 成都 610066；3. 空軍工程大學 防空反導學院，陜西 西安 710051）

0 引言

近年來，隨著窄帶濾光技術和成像光譜技術的快速發展［1-3］，光譜成像已經成為發展天基預警衛星新型探測手段的重要方向［4-6］。光譜圖像包含探測場景的空間信息和光譜信息，具有“圖譜合一”的特點［7］，提供了豐富詳細的物質類別信息，能夠大大提高目標檢測和分類識別精度，但其龐大的數據量也給數據采集和實時處理帶來了巨大挑戰。面對這一矛盾，基于壓縮感知［8-13］的光譜圖像編碼，為實現數據的高效采集提供了一條可行途徑。現有的壓縮采樣編碼多集中在二維空間域，而光譜維則采用常規的數據壓縮方法，忽略了光譜圖像的三維空間稀疏性，仍存在資源浪費，其算法存在循環嵌套，時間復雜度亦為較高的平方階O(n2)。

針對這一問題，本文提出一種基于張量分解的光譜圖像壓縮感知重構方法，該方法利用光譜圖像數據三維空間稀疏性，建立基于三階張量Tucker 分解的光譜圖像重構模型，設計基于正交匹配（orthogonal matching pursuit，OMP）［14］算法的模型求解方法；并將OMP 算法推廣到三維空間，設計一種以三階張量為字典原子的正交匹配追蹤算法，在三維空間實現光譜圖像數據的壓縮采樣及解碼重構。

1 光譜圖像三維空間稀疏性分析

光譜圖像數據通常有光譜空間、圖像空間和特征空間等3 種表達形式［15］，如圖1 所示。

圖1 光譜圖像數據的表示形式Fig. 1 Representation of spectral image data

圖1 中，光譜空間表達為一條連續光譜曲線，包含各像元的光譜信息，反映的是物質輻射特性隨波段變化的關系，不同物質具有不同的光譜曲線；圖像空間表達能夠反映場景中各類物質的空間位置分布和上下文信息；特征空間中，則用多維矢量表示各樣本的光譜向量，能夠定量反映樣本的光譜輻射特性及其變化規律，不同的物質有著不同的統計分布特性。

作為三階張量，光譜圖像包含二維空間數據和一維光譜數據，具有空-譜切片和空-空切片2 種切片分解方式。利用空間域或譜間稀疏性，可以在二維空間域或一維譜域內實現光譜圖像壓縮采樣［16-19］，但實際上光譜圖像三維數據塊在空譜域三維空間內同樣具備稀疏性。以圖2a）為例，構造3個不同的離散小波變換基，基于Tucker 分解計算光譜圖像的變換域核心張量，其向量形式如圖2b）所示。選擇其中10%的顯著系數重構光譜圖像，如圖2c）所示，各波段峰值信噪比（peak signal to noise ratio，PSNR）曲線如圖2d）所示。可以看出，光譜圖像數據在三維空間內具有較強的稀疏性，稀疏分解后，利用部分顯著系數能夠精確重構原始圖像，重構數據各波段的峰值信噪比達到36 dB 左右。可以發現，稀疏分解重構損失了光譜圖像的部分高頻細節信息，導致重構圖像較原始圖像更加平滑，不同類型面目標之間的邊界有一定程度的模糊。

圖2 光譜圖像三維空間稀疏性Fig. 2 Sparsity of spectral image in three-dimensional space

2 基于三階張量分析的壓縮感知光譜圖像重構

基于光譜圖像三維空間稀疏性，直接在三維空間內進行壓縮采樣和解碼重構，可以有效降低計算復雜度，其采樣過程如圖3 所示。場景輻射通過前光學器件后入射到數字微鏡器件（digital mirror device，DMD），空間信息在DMD 上形成一面陣，空間編碼圖像經過聚光透鏡入射到分光計。隨后，分光后的空間編碼數據入射到隨機編碼光電傳感器陣列，進行光譜維采樣編碼，最終得到三維采樣數據。最后，利用RNG 控制開關1-L的開合進行K次0-1 隨機編碼，實現光譜維壓縮采樣編碼，此時其光譜維采樣率為K/L。

圖3 光譜圖像三維空間壓縮采樣Fig. 3 Compressive sensing in three-dimensional space

2.1 TDB-OMP 算法

式中：×n為張量的n-模式積；為變換域系數，即Tucker 分解的核心張量。利用壓縮采樣矩陣和，分別在光譜圖像三維數據的行、列和光譜維進行獨立觀測：

根據Tucker 分解與Kronecker 積的關系，利用向量化算子vec( · )，將式（2）轉化為一維向量形式為

原子更新后，通過求解下列最小二乘問題估計系數向量為

式中：⊙為Khatri-Rao 積［20］；和分別包含了t= 1，2，…，k次循環中從Φ1，Φ2和Φ3中選擇的原子，即vec(Y)為向量化后的壓縮采樣數據。此時，式（6）的解為

式中：( · )#為矩陣的Moore-Penrose 偽逆。式（7）的計算量特別大，利用Khatri-Rao 積性質可以得到

式中：(n= 1，2，3)為第k次循環中新增入Dn的原子。則第k次循環中Z-1的迭代計算式如下：

最終第k次循環結束時，根據式（8）得到系數向量的計算式為

得到系數向量之后，殘差向量為r=y-，將其矩陣化可以得到殘差矩陣：

式中：matr( · )為矩陣化算子。

當滿足中止判據，則中止迭代，否則轉入下一次迭代。經過循環迭代得到變換域系數張量后，利用式（12）可以重構光譜圖像數據：

綜上所述，TDB-OMP 算法通過n-模式積計算更新原子，采用迭代更新的方式計算Zk，有效降低了計算復雜度，其流程如下，其時間復雜度為線性階O(n)。

TDB-OMP 算法

輸入：感知矩陣{Φ1，Φ2，Φ3}，觀測值Y，截斷誤差ε0，稀疏度K0

初始化：

循環：

步驟 1： 計算殘差與感知矩陣的n-模式積，搜索最匹配的原子

步驟 2： 更新標簽集和入選原子集

步驟 3：求解最小二乘問題，估計稀疏系數

步驟 4： 更新殘差

步驟 5： 更新迭代索引k=k+ 1。

End

2.2 3D-OMP 算法

TDB-OMP 算法利用三階張量與感知矩陣的Tucker 積搜索最佳匹配原子，并采用迭代遞推的方式估計稀疏系數，有效降低了算法的計算復雜度。其不足是系數估計和殘差更新仍以向量形式進行，字典中的原子仍是一維向量。為保留光譜圖像數據的三維結構信息，進一步降低計算復雜度和存儲量需求，本節直接利用感知矩陣各原子的外積作為字典原子，將傳統OMP 算法推廣到三維空間，即3DOMP，基于光譜圖像三維空間稀疏性，在三維空間實現光譜圖像壓縮采樣及解碼重構。

式中：wi，j，k為系數三階張量W的元素，定義k=1，2，…，L)為張量原子，其實質是向量和的外積，即。由此，將觀測值Y看成張量原子的線性組合，系數張量W的元素為對應原子的權重系數。觀測值Y到字典中各原子的投影為

由式（14），（15）可知，在3D-OMP 算法中，字典包含N2L個原子，每個原子是維數為M×M×K的三階張量。類似于傳統OMP 算法，每次迭代中3D-OMP 首先從字典中搜索最匹配的原子，并將其添加到已選原子集，然后求解所有原子集對應的權重系數。在第t次迭代中，觀測值Y為已選t個原子的近似線性組合，此時殘差R定義為

權重系數更新的目標是搜索系數向量w=，使得殘差R的l2范數最小，其實質為求解下列最優化問題：

式中：R：：k，Y：：k和分別為三階張量R，Y和Γu正面切片矩陣，則式（17）等效為

式中：tr( · )為矩陣的跡，滿足，因此有

可得

矩陣H的元素的計算如下：

令Y=a?b?c，則的計算式如下：

得到變換域系數張量后，利用式（1）重構原始光譜圖像。綜上，3D-OMP 時間復雜度為線性階O(n)。

其算法如下：

任務：

輸入：感知矩陣{Φ1，Φ2，Φ3}，觀測值Y，截斷誤差ε0，稀疏度T0

初始化：

循環：

步驟 1： 搜索最匹配原子

步驟 3： 分別利用式（21）和式（22）計算H和f

步驟 4： 求解最優化問題式（19），估計稀疏系數

步驟 6：t=t+ 1

End

3 仿真校驗

為檢驗TDB-OMP 和3D-OMP 2 種解碼重構算法的有效性，將它們應用到光譜圖像Our Data［22］，考察它們在不同數據壓縮率下的重構性能，并與文獻［5］中的2D-CRPG+KL 方法進行比較。式（1）中的稀疏基分別為DWT（discrete wavelet transform）變換矩陣Ψ1=Ψ2和DCT（discrete cosine transform）變換矩陣Ψ3。雖然TDB-OMP 和3D-OMP 可以在行、列和光譜維采用任意的測量矩陣組合來保證數據壓縮率，但為確保重構結果的可比性，采用與2D-CRPG 算法相同的測量矩陣組合。給定譜間壓縮率為50%，即K= 150，。假設空間采樣率為Rspa，則空間域采樣矩陣為和，其中。

圖4 給出了空間采樣率為40%～90%時，運用2D-CRPG+KL，TDB-OMP 和3D-OMP 算法重構光譜圖像的曲線。從圖4 可以看出，在相同的數據壓縮率下，TDB-OMP 的重構性能在多數情況下優于2DCRPG+KL，而3D-OMP 算法的重構性能總是優于其他2 種算法，特別是在信噪比較低的波段（波段150～300），其優勢更加明顯。這主要是因為，2DCRPG+KL 對各波段圖像進行獨立壓縮采樣和解碼重構，忽略了波段間的相關信息。TDB-OMP 雖然在三維空間實現光譜圖像重構，使得重構性能有所提升，但其系數和殘差更新仍以向量形式進行，忽略了光譜圖像三維數據的結構信息。盡管遞推求解方式在一定程度上減小了計算量，但TDB-OMP 需要的存儲量仍然較大。而3D-OMP 不僅在三維空間重構光譜圖像，而且直接以三階張量作為字典原子，有效保留了光譜圖像的結構信息，進一步提升了重構性能。同時，由于2D-CRPG+KL 方法在進行譜域KL 變換編碼時，無論譜域壓縮率大還是小都需要計算空間采樣數據的協方差矩陣，并進行特征值分解，使得其計算量和需要的存儲空間都較大。為保證重構質量，3D-OMP 能夠靈活地設置空間采樣率和譜域壓縮率對光譜圖像三維數據進行觀測，而對于2D-CRPG+KL 來說，則需要分析空間采樣率和譜域壓縮率的關系，以及它們對重構性能的影響來優化整體的數據采樣率，如何調整空間采樣率和譜間壓縮率仍然有待解決。

圖4 不同采樣率下重構光譜圖像各波段PSNRFig. 4 PSNRs of reconstructed spectral image with different sampling rate

圖5 給出了2D-CRPG+KL，TDB-OMP，3D-OMP 3 種算法重構光譜圖像的均峰值信噪比（average PSNR，APSNR）與空間采樣率的關系。可以看出，當空間采樣率較小時，TDB-OMP 和3D-OMP 的重構性能明顯優于2D-CRPG+KL，而隨著空間采樣率的增長，3 種算法的APSNR 都逐漸增大。其中TDBOMP 和3D-OMP 都呈近似線性增長，但前者的增長速度小于后者，兩者之間的差距不斷擴大；而2DCRPG+KL 的增長速度較快，當空間采樣率較大時，其重構性能超過TDB-OMP，直至趕上3D-OMP。

圖5 APSNR 隨采樣率變化Fig. 5 Change of APSNR over sampling rate

圖6，7 分別給出了空間采樣率較小（40%）和較大（90%）時，運用3 種算法重構光譜圖像的假色圖。從圖6 可以看出，當空間采樣率較小時，2D-CRPG+KL 重構圖像存在明顯的模糊，而TDB-OMP 重構圖像也與原始圖像有比較明顯的區別，3D-OMP 能夠更好地保留原始圖像的結構信息，其APSNR 略高于TDB-OMP。從圖7 可以看出，當空間采樣率為90%時，不同算法重構光譜圖像的APSNR 均超過39 dB，且2D-CRPG+KL 的APSNR 值最大。這說明，當采樣資源充足，空間采樣率接近全采樣時，3種重構算法均能較好地保留圖像的結構信息和細節特征，三維空間壓縮采樣的優勢并不明顯。

圖6 空間采樣率為40%時不同方法重構光譜圖像Fig. 6 Reconstructed spectral images with 0.4 sampling rate

圖7 空間采樣率為90%時不同方法重構光譜圖像Fig. 7 Reconstructed spectral images with 0.9 sampling rate

圖8 ，9 給出了空間采樣率為40%和90%時，不同算法得到的重構光譜曲線，選取Our Data 中比較典型的3 類樣本，其空間位置分別為（80，45），（170，70）和（151，151）。從圖9 可以看出，當空間采樣率達到90%時，3 種算法均能比較完整地保留典型樣本的光譜信息。而從圖8 可以看出，當采樣率較低時，相比于2D-CRPG+KL 方法，TDB-OMP 和3D-OMP 算法能夠更好地保留各類典型樣本的光譜信息，以更低的資源消耗提升光譜探測信息的應用精度。

圖8 空間采樣率為40%時3 種方法的重構光譜Fig. 8 Reconstructed spectrum with 0.4 sampling rate

圖9 空間采樣率為90%時3 種方法的重構光譜Fig. 9 Reconstructed spectrum with 0.9 sampling rate

4 結束語

本文針對現有光譜圖像壓縮采樣方法仍然存在的高資源消耗問題，利用光譜圖像的三維空間稀疏性和數據塊結構信息，在三維空間建立基于壓縮感知理論的光譜圖像壓縮采樣及重構模型，提出兩種基于三階張量分解的壓縮感知光譜圖像重構算法，有效降低了重構算法計算復雜度，改善了光譜圖像的重構效果，增強了空間采樣率和譜間壓縮率設置的靈活性。