趙梓淵 ，唐意東 ，黃樹彩

（1. 空軍工程大學(xué) 信息與導(dǎo)航學(xué)院，陜西 西安 710077；2. 中國(guó)人民解放軍95607 部隊(duì)，四川 成都 610066；3. 空軍工程大學(xué) 防空反導(dǎo)學(xué)院，陜西 西安 710051）

0 引言

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

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

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

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

圖1 光譜圖像數(shù)據(jù)的表示形式Fig. 1 Representation of spectral image data

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

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

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

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

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

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

2.1 TDB-OMP 算法

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

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

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

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

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

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

最終第k次循環(huán)結(jié)束時(shí)，根據(jù)式（8）得到系數(shù)向量的計(jì)算式為

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

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

當(dāng)滿足中止判據(jù)，則中止迭代，否則轉(zhuǎn)入下一次迭代。經(jīng)過循環(huán)迭代得到變換域系數(shù)張量后，利用式（12）可以重構(gòu)光譜圖像數(shù)據(jù)：

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

TDB-OMP 算法

輸入：感知矩陣{Φ1，Φ2，Φ3}，觀測(cè)值Y，截?cái)嗾`差ε0，稀疏度K0

初始化：

循環(huán)：

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

步驟 2： 更新標(biāo)簽集和入選原子集

步驟 3：求解最小二乘問題，估計(jì)稀疏系數(shù)

步驟 4： 更新殘差

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

End

2.2 3D-OMP 算法

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

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

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

權(quán)重系數(shù)更新的目標(biāo)是搜索系數(shù)向量w=，使得殘差R的l2范數(shù)最小，其實(shí)質(zhì)為求解下列最優(yōu)化問題：

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

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

可得

矩陣H的元素的計(jì)算如下：

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

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

其算法如下：

任務(wù)：

輸入：感知矩陣{Φ1，Φ2，Φ3}，觀測(cè)值Y，截?cái)嗾`差ε0，稀疏度T0

初始化：

循環(huán)：

步驟 1： 搜索最匹配原子

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

步驟 4： 求解最優(yōu)化問題式（19），估計(jì)稀疏系數(shù)

步驟 6：t=t+ 1

End

3 仿真校驗(yàn)

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

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

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

圖5 給出了2D-CRPG+KL，TDB-OMP，3D-OMP 3 種算法重構(gòu)光譜圖像的均峰值信噪比（average PSNR，APSNR）與空間采樣率的關(guān)系?？梢钥闯?，當(dāng)空間采樣率較小時(shí)，TDB-OMP 和3D-OMP 的重構(gòu)性能明顯優(yōu)于2D-CRPG+KL，而隨著空間采樣率的增長(zhǎng)，3 種算法的APSNR 都逐漸增大。其中TDBOMP 和3D-OMP 都呈近似線性增長(zhǎng)，但前者的增長(zhǎng)速度小于后者，兩者之間的差距不斷擴(kuò)大；而2DCRPG+KL 的增長(zhǎng)速度較快，當(dāng)空間采樣率較大時(shí)，其重構(gòu)性能超過TDB-OMP，直至趕上3D-OMP。

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

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

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

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

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

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

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

4 結(jié)束語

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