面向誤差補償的高光譜與多光譜圖像融合

方帥，許漫

1.合肥工業大學計算機與信息學院,合肥 230601；2.工業安全與應急技術安徽省重點實驗室，合肥 230601

0 引 言

高光譜圖像(hyperspectral image, HSI)提供豐富的光譜信息，廣泛應用于圖像分類和目標檢測等領域。然而，現有成像的光譜分辨率和空間分辨率是相互制約的。因此，HSI通常具有較低的空間分辨率，以保證較高的信噪比，然而如此低的空間分辨率限制了HSI的適用性。與高光譜圖像相比，多光譜圖像(multispectral image, MSI)通常具有更高的空間分辨率，但是光譜波段較少。將同一場景的低空間分辨率的高光譜圖像與高空間分辨率的多光譜圖像進行融合，可以重建高空間分辨率的高光譜圖像(high-spatial-resolution hyperspectral image, HR-HSI)。

大多現有的方法受光譜字典學習算法啟發，先從HSI中學習光譜字典，后對MSI進行稀疏編碼來估計系數，從而獲得重建的高分辨率高光譜圖像。然而僅僅利用光譜字典進行重建目標圖像會不可避免地損失一部分空間信息。同時由于遙感圖像的地物分布極為復雜，因此其空間信息在不同的區域可能是千差萬別的?，F有的大多進行誤差補償的算法(Yang等，2020；Yin和Li，2015)也僅考慮了不同通道之間的區別，而忽略了空間信息分布的復雜性，因而存在局限性。

圖1展示了Pavia University數據的波段1、波段30、波段50和波段70的水平方向梯度及局部放大圖，第1行分別為4個波段的梯度，紅色框中區域的放大圖對應第2行，藍色框中區域的放大圖對應第3行。從圖1中可以看出，不同通道在同一區域表現出了不同的結構特性，同一通道在不同區域也表現出了不同的結構特性。因此，進行誤差注入時若只考慮通道之間的不同則很容易忽略每個通道的不同區域之間的差異，從而導致融合結果質量降低。

圖1 Pavia University數據梯度及其局部放大圖

本文提出了通過一種基于誤差補償的算法來重建高分辨率高光譜圖像，即先通過光譜字典學習算法來獲得中間結果，再利用殘余空間信息對中間結果進行補償，使得融合結果的質量進一步提高，本文主要貢獻總結如下：

1)針對光譜字典學習算法損失空間信息這一問題，提出了誤差補償策略，先利用低分辨率圖像構建出空間殘差，后在光譜字典學習的基礎上進行空間殘差信息注入，設計了一種局部區域注入系數算法，使得到的系數能夠根據局部區域的光譜特征自適應地調節誤差信息的注入，從而避免注入過多或過少的空間信息而導致光譜失真。

2)為保證經過誤差補償后的融合結果在提升空間分辨率的同時不發生光譜畸變，在梯度域探索MSI和HR-HSI之間的關系，構建變分模型。并將該變分模型與1)中構成的優化項結合，從而構建出一個新的模型，該模型中既含有空間方面的約束又含有光譜方面的約束，通過迭代求解逐步提升系數的精度和融合圖像的質量。

1 相關工作

盡管各個應用領域對于高分辨率高光譜圖像有廣泛需求，但由于技術和成本的限制，仍難以同時實現空間高分辨率和光譜高分辨率，因此這個問題引起了廣泛關注。學者們利用現有的圖像融合技術將同一場景下的HSI和MSI/RGB(red-green-blue)進行融合，從而重建出HR-HSI。

早期高光譜與多光譜圖像(HSI-MSI)融合算法從Pansharpening(Alparone等，2007；方帥 等，2020；焦姣和吳玲達，2019)算法借鑒而來，該類算法側重于細節的提取和注入，由于全色圖像不包含任何光譜信息，因此，融合結果通常存在嚴重光譜畸變。隨后出現了一批獨立發展起來的HSI-MSI融合算法，主要可分為3類：基于貝葉斯的方法、基于深度學習的方法和基于低維子空間的方法。

基于貝葉斯的方法(Hardie等，2004；Akhtar等，2015)利用觀測場景中的先驗分布，如高斯先驗，并結合最大后驗概率實現對問題進行求解。Wei等人(2015)提出一種快速融合的方法，其通過求解與貝葉斯估計相關的Sylvester方程達到融合目的。但是貝葉斯算法依賴于先驗假設，而人工設計的先驗假設通常是針對某一方面的考慮，并不全面，因而可能導致融合圖像質量的降低。

由于深度卷積神經網絡(deep convolutional neural networks, DCNN)(Dong等，2014)能夠對提取的特征進行分析學習，在HSI-MSI融合領域展現出優越的性能，從而得到廣泛應用。Yang等人(2018)提出兩分支卷積神經網絡(convolutional neural networks,CNN)結構，分別從輸入的HSI和MSI中提取光譜特征和空間特征，然后將特征拼接并傳遞給全連接層，生成HR-HSI。Dian等人(2018b)提出了一種基于深度殘差學習的HSI超分辨率方法，該方法使用DCNN學習初始化的HR-HSI和ground truth之間的映射關系。但是基于深度學習的算法通常需要大量的數據用于訓練模型，而當前高光譜圖像來源較少，這一問題限制了深度學習在高光譜圖像融合中的應用。Xie等人(2022)提出了一種MHF-net(MS/HS fusion net)，該方法使用DCNN去噪器估計HR-HSI中間結果的正則化，從而進一步用于求解具有低秩約束的目標HR-HSI。

基于低維子空間的方法假設低分辨率高光譜圖像或高分辨率高光譜圖像的光譜像素存在于低維子空間，通過提取光譜特征和求取系數來重建HR-HSI，這類方法可分為基于張量分解的方法和基于字典學習的方法?；趶埩糠纸獾姆椒ㄓ靡粋€核心張量和3個模方向的字典乘積表示目標高空間分辨率高光譜圖像。3D張量可以很好地表示同樣具有3維結構的HSI和MSI，基于此，張量分解廣泛應用于HSI-MSI融合問題中。Xu等人(2022)提出一種基于張量子空間表示的迭代正則化方法來進行HSI-MSI融合，該方法不僅應用了HR-HSI的全局譜空間低秩與非局部自相似性先驗，而且在迭代過程中利用融合結果與輸入的低分辨率圖像之間的殘差信息進行誤差補償。但是該方法在進行迭代補償殘差信息時僅將殘差做同等處理，而并未考慮空間信息的補償，因此方法性能還有很大的提升空間。與矩陣的秩不同，張量的秩是一個非確定性多項式(non-deterministic polynomial, NP)問題,因此基于張量分解的方法通常面臨秩的不確定性和分解的不唯一性等問題。

基于字典學習的方法包括空間字典學習和光譜字典學習，對于HSI-MSI融合問題來說，光譜信息比空間信息重要的多，因此光譜字典學習的使用更為廣泛。Kawakami等人(2011)首先將稀疏矩陣分解(sparse matrix factorization, SMF)用于HSI-MSI融合，先用稀疏字典學習法估計出光譜字典，后從MSI中估計出系數。Akhtar等人(2014)和Dong等人(2016)設計了一種利用HSI的空間信息非負性來約束解空間的稀疏表示(sparse representation, SR)方法。Han等人(2018)和Fang等人(2018)提出了基于超像素的SR模型，該模型根據空間局部結構自適應調整聚類簇的大小和形狀，提高了算法的精度和效率。Dian等人(2019)將非局部相似性與稀疏性先驗相結合，同樣提高了算法的精度。這幾種方法是基于字典學習算法中較為經典的方法，利用先驗信息構建正則項來盡可能捕獲空間信息，然而這并不足以充分捕獲HSI和MSI中的空間信息。為了解決這一問題，Han等人(2020)提出了基于雙字典的高光譜圖像的融合方法，雙字典包括一個常規的光譜字典和一個針對殘差圖的空間字典，空間字典是為了表征光譜字典重構融合圖像而導致空間信息的損失，該方法由于兼顧了光譜信息和空間信息的表示，具有優秀的融合性能。

本文方法與光譜字典學習算法有關，但又與其有很大差異。本文只將字典學習算法所得結果作為中間結果，通過退化模型構造誤差對其進行補償，使得融合結果質量進一步提升。在進行空間殘余信息注入時，設計了一種局部區域注入系數，使得到的系數能夠根據局部區域的光譜特征自適應地調節誤差信息的注入。在殘差信息注入的同時還考慮了保持光譜信息不發生畸變，通過優化目標函數，使得到的融合結果在空間和光譜性能上均表現出一定的優越性。

2 問題描述

高光譜圖像和多光譜圖像融合算法可描述為通過融合同一場景下的HSI和MSI獲得HR-HSI。HSI、MSI以及HR-HSI分別用X∈RL×n、Y∈Rl×N和Z∈RL×N表示，其中L和l(L>l)表示圖像的波段數，N和n(N>n)表示圖像的像素點數目,則X和Y可以表示為Z的線性變換，即

X=ZBM+E1

(1)

Y=FZ+E2

(2)

式中，矩陣B∈RN×N是點擴散函數，M∈RN×n是一個下采樣矩陣,F∈Rl×L表示Z映射到Y的光譜響應函數。E1、E2表示殘差。

如果用稀疏表示對圖像Z進行描述，可表示為

Z=DA+E=T+E

(3)

式中，D∈RL×K表示光譜字典；A∈RK×N是對應的稀疏系數矩陣；K表示字典的原子數，N表示圖像像素點數，且K?N；E∈RL×N表示不能被光譜字典D表達的誤差矩陣。

由式(3),結合式(1)和式(2)，有

X=DABM+EH

(4)

Y=FDA+EM

(5)

式中，EH∈RL×n,EM∈Rl×N分別代表HSI和MSI中的誤差矩陣，其表示不能被光譜字典表達的空間信息，可以從誤差矩陣EH和EM中挖掘出需要的空間殘差信息。由于EM中含有更多的空間信息，因此后續主要對EM進行處理。

在進行誤差補償時，首先將初步結果圖像T和多光譜圖像Y按照相關性進行分組，即計算T的所有波段與Y每個波段的相關系數，得到一個相關系數表。根據表中每列的最大值將T當前波段劃分到對應的Y的一組內，對T的所有波段都照此進行，Y有多少個波段就將T分為多少組。將每一組的融合結果Zj(j=1,2,…,l)表示為

Zj=Tj+g.*(Yj-I)

(6)

式中，Tj(j=1,2,…,l)表示第j組圖像T；Yj(j=1,2,…,l)表示多光譜圖像的第j個波段；g表示誤差補償系數；“.*”表示點乘運算；I表示對于圖像T模擬生成的低分辨率亮度分量，這里對圖像T進行光譜下采樣直接得到I，有

Zj=Tj+g.*(Yj-(FT)j)=

(7)

本文的目標是精確地求解光譜字典D和系數矩陣A以得到圖像T，以及在保證不發生光譜畸變的情況下求得誤差矩陣補償系數g進行誤差補償，從而獲得精確穩定的融合結果。

3 融合算法

本文主要包括兩個模塊：初步結果估計模塊以及誤差補償模塊,如圖2所示。初步結果估計模塊包括光譜字典訓練和基于超像素分割的系數估計。通過輸入圖像的局部相似性構造正則項，縮小解空間以估計出稀疏系數，從而獲取精確穩定的解。誤差補償模塊包括誤差補償系數的估計和目標圖像Z的重建。設計了一種局部區域注入系數，使得到的補償系數更能反映對應場景的特征，且為了保證經過誤差補償后的融合結果在提升空間分辨率的同時不發生光譜畸變，在梯度域探索融合圖像Z和多光譜圖像MSI的關系，構建一個新的目標函數，設計了一種迭代方法，以逐步提高注入系數的精度和目標圖像Z的質量。

圖2 本文算法框架圖

3.1 初步結果估計

3.1.1 字典訓練

由于融合圖像Z的光譜信息與HSI的光譜信息一致，因此從HSI中學習光譜字典D。由于HSI中具有大量光譜特征相似的地物目標，依據光譜特性采用K-means對HSI進行聚類，再針對每類學習光譜字典，保障字典對各類目標都具有可靠的表達能力，有

i=1,2,…,p

(8)

3.1.2 系數估計

一個典型的自然場景通常包含圖像中所有相似像素的集合，而在一個鄰域內的像素可能表示相同的地物，因此HR-HSI是局部低秩的。當局部區域變化不豐富時滿足局部低秩，為了使用局部低秩先驗，首先使用SLIC(simple linear iterative clustering)(Wang等，2017)超像素分割算法對MSI進行分割，在超像素塊空間內像素具有同質特性，從而滿足低秩假設。同理，HR-HSI也被分成同等數目的超像素，然后在目標函數中添加低秩正則項，構成一個新的目標函數。

局部相似約束可以顯著抑制噪聲，從而獲得更可靠的融合結果。為了利用HR-HSI的局部相似性，假設一個像素可以由與其相似的像素的線性組合來近似，換句話說，該像素的稀疏向量可由同一超像素內所有像素的稀疏向量加權平均得到，即

(9)

構建以下目標函數進行系數求解

(10)

式(10)是凸函數，因此可以用交替方向乘子法(alternative direction multiplier method, ADMM)(Boyd等，2011)進行有效求解。

3.2 誤差補償

3.2.1 誤差補償系數估計

由問題描述中的分析可知，本文任務之一就是如何將誤差EM精確地補償到初步結果T中。設計一種局部區域誤差注入的模型，即利用超像素分割算法將圖像分割成多個區域，同一超像素塊中的像素具有同質特性，相同的系數增益應用于同一超像素塊中的所有像素，從而能夠根據局部區域的光譜特性自適應地調整誤差信息的注入權重。該模型是一種局部光譜驅動的注入模型，能夠盡可能地保持融合后的目標圖像的光譜信息不產生畸變。在進行誤差補償時，將初步結果圖像T和多光譜圖像Y進行分組，每一組對應的補償系數為

(11)

現有的常見方法都只是針對不同的波段計算不同的g，這些方法忽略了圖像空間信息的分布是很復雜的，同一波段的不同區域中的結構特性也會相差很大，從圖1中即可看出。因此細節恢復并不是只關注通道變化這么簡單,同時也證明了本文設計的系數合理性。

3.2.2 目標圖像Z估計

現有誤差注入算法在進行誤差注入時只關注提升空間分辨率，而幾乎不會考慮注入過程中光譜信息是否發生畸變，因此融合結果的性能還有很大的提升空間。本文不是僅僅依靠式(7)來獲得最終融合結果，而是進一步通過構造一個同時包含空間信息約束和光譜信息約束的目標函數來獲得最終的融合結果，使得進行誤差補償后的融合結果在提升空間分辨率的同時不發生光譜畸變。構造以下目標函數來實現，即

(12)

(13)

式中，F為光譜響應函數。

式(12)的3個優化項都是凸函數，因此同樣可以使用ADMM對其進行有效求解。

為了提升補償系數g的精度和融合結果Z的質量，通過迭代更新式(11)和式(12)來實現。即先初始化g為

(14)

式中，T為字典學習算法所得初步結果，然后用初始化的g通過式(12)計算出Z，再用計算出的Z重新計算g，即

(15)

重復上述迭代過程，直至達到設定的迭代次數。通過這樣g與Z相互迭代更新，由于融合結果Z的質量越來越高，使得迭代過程中g的估計越來越精確。又正是因為g的估計越來精確，使得補償的誤差越來越精準，從而使得Z的質量越來越高，這是一個相互促進的過程。

4 實驗與分析

4.1 實驗數據與參數設置

4.1.1 Pavia University數據集

4.1.2 AVIRIS數據集

第2個數據集是由機載可見光/紅外成像光譜儀(airborne visible infrared imaging spectrometer, AVIRIS)獲取，圖像尺寸為300×300像素。數據集的波段1—2，105—115，150—170，223—224等被水蒸氣吸收波段后，剩下93個波段，類IKONOS的反射光譜響應濾波器用于生成HR-MSI。使用模板為5×5、標準差為3的高斯濾波器對ground truth進行高斯濾波，然后對每個波段在水平方向和垂直方向上分別進行4倍下采樣得到LR-HSI。

4.1.3 真實數據集

4.1.4 參數設置

圖3 隨參數變化的PSNR曲線

4.2 對比算法和評估標準

為了驗證算法的整體有效性，選擇6種優秀算法與本文算法進行對比，分別是CNMF(coupled nonnegative matrix factorization)(Yokoya等，2012)、Hysure(Sim?es等，2015)、NSSR(non-negative structured sparse representation)(Dong等，2016)、LRSR(low-rank and sparse representations)(Dian等，2018a)、OTD(optimized twin dictionaries)(Han等，2020)和IR-TenSR(iterative regularization method based on tensor subspace representation)(Xu等，2022)，其中NSSR、LRSR和OTD是基于字典學習的融合算法。

為了對各算法融合結果進行評價，本文選取4種評價指標檢測融合結果的光譜質量和空間質量，分別是峰值信噪比(PSNR)、全局相對誤差(relative dimensionless global error in synthesis, ERGAS)、光譜角制圖(spectral angle mapper, SAM)和均方根誤差(root mean square error, RMSE)。

4.3 實驗結果與分析

4.3.1 實驗設計

為了驗證本文基于誤差補償的高光譜與多光譜圖像融合算法的有效性，設計以下實驗：

實驗1)：為了驗證本文誤差補償模塊的有效性，在其他優化項以及正則化參數一致的情況下，分別建立初步結果估計模塊結果和誤差補償后的結果，對比使用誤差補償前后對融合結果的影響。

實驗2)：為了驗證本文局部區域誤差補償系數估計的有效性，采用相同方式獲得字典學習結果，且在誤差補償模塊保證其他優化項相同的基礎上，對比局部注入系數與全局注入系數對融合結果的影響。

實驗3)：為了驗證本文算法有效性，在2個模擬數據集以及1個真實數據集上與其他6個算法進行比較。

實驗4)：為了進一步驗證本文算法誤差補償部分的普適性，對各對比算法增加誤差補償作為后處理，對比誤差補償算法對于融合效果的影響。

所有算法均在英特爾酷睿i5-8400 CPU和8 GB RAM且裝有MATLAB2018a計算機上運行。兩個仿真數據集上采用與IKONOS衛星相同的光譜響應函數F(Dian等，2018a)。此外，光譜字典原子數均設置為100。

4.3.2 實驗結果與分析

實驗1)對比有無誤差補償對融合效果的影響，圖4分別為加上誤差補償模塊前和加上誤差補償模塊后重構的HR-HSI與ground truth之間的誤差圖及其局部區域放大圖。對比可以發現，加上誤差補償后融合結果產生的畸變更小，圖像的質量有所提高。從表1中可以看出，進行誤差補償后，融合結果的各項指標均優于無誤差補償的結果，尤其在光譜保持方面有很明顯的效果。圖5為各算法在AVIRIS數據集上的SAM圖，從圖5中可以看出,本文算法所得結果的光譜畸變最小，由此說明本文方法的融合結果的光譜特性最佳。

圖4 誤差圖

表1 誤差補償項對比實驗結果

實驗2)在均使用誤差補償模塊的情況下，對比局部注入系數與全局注入系數對融合效果的影響。從表2可以看出，局部注入系數無論是在空間結構方面還是光譜信息保持方面都優于全局注入系數，因此利用局部注入系數所得融合結果的性能更為優越。

表2 注入系數對比實驗結果

實驗3)在與各算法的對比實驗中，圖6和圖7分別展示了Pavia University數據集第10、16、47波段的偽彩色融合結果和AVIRIS數據集第49、24、11波段的偽彩色融合結果及其對應的MSE圖。第1行為融合結果，第2行為對應的MSE圖，第3行為MSE圖的局部放大圖，放大區域為MSE圖中紅色框部分。從上述圖中可以看出，盡管所有算法的融合結果與ground truth相比都存在一些畸變，但是都能恢復大部分空間細節。其中與CNMF、Hysure、NSSR和LRSR這些只利用光譜基或光譜字典的算法比較而言，由于IR-TenSR利用殘差信息進行誤差補償，充分利用低分辨率圖像的信息，因此誤差較??；而本文算法與OTD算法不僅僅使用光譜字典來重建HR-HSI，還利用殘差空間信息進行誤差補償，因此本文算法和OTD算法保持空間結構的性能也很優越。

圖6 Pavia University數據實驗結果及對應MSE圖

圖7 AVIRIS數據實驗結果及對應MSE圖

圖8為各算法在真實數據集上的融合結果，由于真實數據集無參考圖像，因此給出其在第16、5、2波段的偽彩色圖像，其中圖8(a)為LR-HSI第16、5、2波段的偽彩圖，圖8(b)為HR-MSI第1、2、3波段的偽彩圖，由于CNMF融合結果過差，因此未展示CNMF融合結果。從圖8中可以看出，除CNMF外，其他算法均能提高LR-HSI的分辨率，但與HR-MSI相比，Hysure和NSSR，LRSR算法不能很好地恢復空間結構, OTD算法融合后的圖像存在嚴重的光譜失真，而IR-TenSR算法光譜失真較小但是空間分辨率較低。相比之下，盡管本文算法重建圖像不完美，但是其在光譜失真和空間結構保留之間取得了平衡。

圖8 真實數據集實驗結果

表3和表4中各指標名稱下左欄為各算法融合結果在兩個數據集上的質量評價標準下獲得的指標值。從客觀的指標值中不難看出，由于本文算法在提升空間分辨率的同時將保持光譜信息不發生畸變考慮在內，因此本文算法在空間和光譜性能上均表現出一定的優越性。由于OTD算法在注入空間誤差時未考慮空間信息分布的復雜性而直接進行全局注入，導致其融合結果依然存在提升空間。從表3中可以看出，本文算法相對于OTD算法在SAM指標上的精度分別提升了3.5%和8.1%，這對于高光譜圖像來說非常重要，因此可以看出本文在光譜信息保持方面有較好效果。

實驗4)的結果可參見表3和表4。表3和表4中各指標名稱下后處理一欄為各對比算法加上本文算法誤差補償部分后的定量評價指標值。比較表3和表4中兩欄指標結果可以看出，在加上誤差補償部分作為后處理后，各算法融合結果的質量都有一定的提升，這也說明了本文誤差補償策略的有效性與價值意義。

表3 Pavia University數據實驗結果

表4 AVIRIS數據實驗結果

從表5可以看出，本文算法在時間性能上沒有優勢，這是因為本文在字典學習算法的基礎上添加了誤差補償策略，具有兩個目標函數需要進行優化更新，這使得本文的時間復雜度增加、運行時間增長。OTD算法與本文算法類似，均進行空間誤差補償，其算法運行時間與本文相當。

表5 測試算法運行時間

5 結 論

通過對現有高光譜與多光譜圖像融合算法的分析，本文提出一種基于誤差補償的融合算法，利用高光譜與多光譜之間的成像關系，計算重建誤差在多光譜的各個通道的體現；將多光譜誤差以不同比例對高光譜的各通道和各個區域進行精準補償。在進行誤差補償時，提出基于局部區域的誤差補償系數，使得到的補償系數更能反映對應場景的特征，因此能夠根據局部區域的光譜特性自適應地調整誤差信息的注入權重。同時，為了保證經過誤差補償后的融合結果在提高分辨率的同時保持光譜信息不發生畸變，在梯度域提取MSI的空間結構，構造變分模型，并將該變分模型與誤差補償系數優化項結合，進行迭代更新。在Pavia University和AVIRIS兩個數據集上，與當前具有代表性的6種優秀融合算法進行比較，本文算法在評價指標以及視覺效果上都取得了有效的提升；在無參考的真實數據集上，本文算法結果在空間和光譜上都呈現出更好的視覺效果。

本文算法主要解決基于光譜字典的融合算法在空間信息誤差較大的問題。利用光譜字典得到初步融合結果，而重點是對空間誤差的補償。假設初步融合結果已經充分獲得光譜信息，在進行誤差補償時盡量保持光譜信息不發生變化。然而,忽略初步結果在光譜方面同樣存在誤差的問題，因此，本文算法在光譜信息方面還有很大的提升空間。在未來的研究中，需要提高光譜信息的準確性。此外還可以在不同的域(而不僅僅是梯度領域)中探索HSI和MSI之間的關系。