蔣彤，安如*，邢菲，王本林, 琚鋒

(1.河海大學(xué) 水文水資源學(xué)院，南京 210098;2.河海大學(xué) 地球科學(xué)與工程學(xué)院，南京 211100;3.滁州學(xué)院 地理信息與旅游學(xué)院，安徽 滁州 239001)

0 引言

高光譜影像擁有較多波段，能夠準(zhǔn)確地反映不同地物特有的吸收或反射特性[1-2]。但由于受到成像機(jī)制的限制，高光譜影像通?？臻g分辨率較低，混合像元分布嚴(yán)重，且幅寬較小，在應(yīng)用中有著局限性[3]。多光譜影像雖然波段數(shù)少，但幅寬較大，且空間細(xì)節(jié)信息豐富。針對(duì)高光譜影像和多光譜影像各自的優(yōu)勢(shì)，通過(guò)圖像融合算法，將多光譜影像中的空間細(xì)節(jié)信息與高光譜影像中的光譜信息進(jìn)行融合，獲得同時(shí)兼具高光譜與高空間分辨率的融合影像，從而提高地物識(shí)別與特征提取能力[4]。

目前，常見(jiàn)的多光譜-高光譜融合算法主要分為3類。一是頻率域?yàn)V波方式。二是基于概率統(tǒng)計(jì)理論的影像融合[5]。前者由于空間信息來(lái)源不一致，導(dǎo)致現(xiàn)有的高光譜-全色的融合方法適用性低[6]。后者在使用梯度下降法求解時(shí)，存在魯棒性不高，計(jì)算復(fù)雜的問(wèn)題[7]。三是最常見(jiàn)的基于混合像元分解理論進(jìn)行數(shù)據(jù)融合。例如Yokaya等提出的聯(lián)合非負(fù)矩陣分解方法(Coupled Nonnegative Matrix Factorization, CNMF)[8]。但現(xiàn)有的像元分解算法只利用了混合像元的光譜維信息，沒(méi)有考慮像元的空間相關(guān)性[9]?；诜秦?fù)矩陣分解的算法由于目標(biāo)函數(shù)非凸性，對(duì)初值敏感，導(dǎo)致端元與豐度矩陣不穩(wěn)定[10]，影響融合圖像的質(zhì)量。

上述三類常見(jiàn)的融合算法大多要求多光譜和高光譜影像完全重疊，無(wú)法對(duì)重疊范圍之外的區(qū)域進(jìn)行融合。Sun等人提出的多光譜圖像光譜分辨率增強(qiáng)方法(Spectral Resolution Enhancement Method, SREM)[10]，實(shí)現(xiàn)了融合方法在空間上的擴(kuò)展。根據(jù)重疊區(qū)高光譜純凈端元，進(jìn)行支持向量機(jī)分類 (Support Vector Machine, SVM)并提取高光譜與多光譜數(shù)據(jù)集，計(jì)算轉(zhuǎn)換矩陣，通過(guò)光譜加權(quán)最小距離篩選轉(zhuǎn)換矩陣，實(shí)現(xiàn)高光譜圖像重構(gòu)，重構(gòu)出的高光譜圖像同時(shí)具有高光譜與高空間分辨率，且增加了幅寬。但在融合過(guò)程中，由于對(duì)圖像進(jìn)行了重采樣和空間配準(zhǔn)，較大的分辨率差、影像不同灰度屬性或?qū)Ρ榷仁沟门錅?zhǔn)出現(xiàn)困難，對(duì)像元的位置和信息量造成影響[11]。

針對(duì)以上問(wèn)題，本研究提出了改進(jìn)的多光譜圖像光譜分辨率增強(qiáng)方法(Improved Spectral Resolution Enhancement Method-Spectral Gradient Angle, ISREM-SGA)，其主要思想是通過(guò)計(jì)算高光譜影像上像元光譜，與多光譜影像上對(duì)應(yīng)位置及其鄰域內(nèi)像元光譜的梯度角，匹配同名點(diǎn)，來(lái)獲得最具光譜相似特性的光譜數(shù)據(jù)集，進(jìn)行高光譜圖像的重構(gòu)。本研究以意大利帕維亞大學(xué)Pavia U(Pavia University)、 日本筑西市Chikusei的高光譜圖像為例，對(duì)兩組原始高光譜數(shù)據(jù)進(jìn)行空間和光譜上的降采樣生成模擬的高光譜、多光譜實(shí)驗(yàn)數(shù)據(jù)進(jìn)行融合研究，最后對(duì)融合圖像進(jìn)行質(zhì)量評(píng)價(jià)。

1 研究數(shù)據(jù)

本研究選取兩組空間變異度差距較大的高光譜圖像為研究數(shù)據(jù)。數(shù)據(jù)集1為2014年7月獲得的Chikusei高光譜圖像，如圖1a(1)所示，由美國(guó)Headwall公司通過(guò)可見(jiàn)光-近紅外C系列高光譜傳感器拍攝生成。此次研究只選用400×400大小的圖像作為基礎(chǔ)數(shù)據(jù)，圖像的空間分辨率為2.5 m，波長(zhǎng)范圍為 0.360～1.018 μm，共128個(gè)波段。圖像上共有部分水體、7種植被、3種裸土以及8種人造物這19類地物[12]，但圖像主要以耕地為主，空間變異度低，光譜特征差異小。Chikusei 數(shù)據(jù)下載地址：http://park.itc.u-tokyo.ac.jp/sal/hyperdata/。數(shù)據(jù)集2選用的是2003年成像的意大利北部帕維亞大學(xué)Pavia U數(shù)據(jù)，如圖1b(1)所示。該高光譜數(shù)據(jù)是由機(jī)載反射光學(xué)光譜成像儀傳感器拍攝生成，共包括從0.43～0.86 μm的103個(gè)波段，每個(gè)波段上像元總數(shù)為610×340，圖像的空間分辨率為1.3 m，圖像包含瀝青道路、磚石、金屬板、柏油房頂、草地等9類地物[12]。數(shù)據(jù)集2是以居民地和交通用地為主，地物類型少，但各光譜特征之間差異較大。Pavia U下載地址：http://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes。

圖 1 實(shí)驗(yàn)數(shù)據(jù)Figure 1 Experimental Data

表 1 Chikusei與Pavia U數(shù)據(jù)多光譜波段信息Table 1 Multispectral band information of Chikusei and Pavia U

對(duì)數(shù)據(jù)集1和2進(jìn)行光譜降維，利用光譜變換矩陣模擬出多光譜圖像，如圖1a(2)和b(2)所示。光譜變換矩陣根據(jù)多光譜圖像各波段的光譜響應(yīng)函數(shù)得出，兩組數(shù)據(jù)均以環(huán)境與災(zāi)害監(jiān)測(cè)預(yù)報(bào)小衛(wèi)星A星(HJ-1A) CCD多光譜數(shù)據(jù)的光譜響應(yīng)函數(shù)為參考，具體操作如下：(1)對(duì)每一條光譜響應(yīng)函數(shù)在原始高光譜數(shù)據(jù)中心波長(zhǎng)處進(jìn)行重采樣，并計(jì)算采樣值在各中心波長(zhǎng)處的權(quán)重；(2) 根據(jù)權(quán)重對(duì)高光譜圖像的波段數(shù)進(jìn)行加權(quán)求和，得到模擬多光譜圖像。模擬出的多光譜圖像擁有4條波段，各波段的詳細(xì)范圍展示在表1中。

對(duì)兩組數(shù)據(jù)集的原始高光譜數(shù)據(jù)進(jìn)行降采樣，獲得空間分辨率較低的模擬高光譜數(shù)據(jù)，對(duì)模擬結(jié)果進(jìn)行裁剪，獲得與多光譜圖像僅有部分重疊的高光譜圖像用于融合實(shí)驗(yàn)。為保證提取到相同位置的像元光譜，利用最近鄰法將高光譜圖像重采樣成多光譜圖像大小，并進(jìn)行空間配準(zhǔn)，結(jié)果展示在上圖1a(3)、b(3)中。但Chikusei 、Pavia U高光譜圖像在前5個(gè)波段出現(xiàn)下圖2中a(1)、b(1)所示的明顯噪聲點(diǎn)，a(3)、b(3)信噪比結(jié)果圖中也表明其信噪比略低，因此只選用剩余無(wú)明顯噪聲點(diǎn)的波段作為最終融合實(shí)驗(yàn)所用的高光譜圖像。

2 研究方法

2.1 多光譜圖像光譜分辨率增強(qiáng)方法:SREM

由于高光譜圖像與多光譜圖像在重訪周期和幅寬上都有著明顯的差異，大多數(shù)融合方法受到覆蓋范圍的限制，但盡管如此，不同傳感器拍攝的圖像上同類地物在光譜信息上仍具有一定的相似性?；谶@一原理， Sun等人提出了SREM多光譜圖像光譜分辨率增強(qiáng)方法。在SREM方法中首先利用頂點(diǎn)成分分析(Vertex Component Analysis, VCA)算法獲得高光譜圖像端元，并根據(jù)純凈端元對(duì)高光譜圖像進(jìn)行SVM分類，根據(jù)分類結(jié)果提取在空間位置上相互對(duì)應(yīng)的各地物高光譜數(shù)據(jù)集、多光譜數(shù)據(jù)集。對(duì)于有著N類地物的圖像而言，可分別用下式(1～2)的矩陣M(q)、H(q) 來(lái)描述第q類地物的多光譜與高光譜數(shù)據(jù)集：

(1)

(2)

式(1～2)中：M(q) 是維度為L(zhǎng)×W的矩陣，H(q) 是維度為K×W的矩陣L與K分別表示多光譜和高光譜波段數(shù)，W表示提取像元的個(gè)數(shù)，W的值不得小于多光譜的波段數(shù),否則無(wú)法提取足量用于影像重構(gòu)的光譜和空間信息。

圖 2 高光譜圖像信噪比與噪聲波段Figure 2 SNR and noisy band of hyperspectral image

G(q)=H(q)M(q)T(M(q)M(q)T)-1

(3)

h(q)=G(q)×M(q)

(4)

其次計(jì)算轉(zhuǎn)換矩陣。Winter等人在提出的CRISP方法[14]中表示，可以通過(guò)移動(dòng)窗口提取光譜數(shù)據(jù)集，并根據(jù)最小二乘方法即公式(3)計(jì)算窗口對(duì)應(yīng)的轉(zhuǎn)換矩陣，將得到的轉(zhuǎn)換矩陣按照公式(4)計(jì)算重構(gòu)高光譜影像。

最后對(duì)轉(zhuǎn)換矩陣進(jìn)行逐像元篩選，實(shí)現(xiàn)高光譜影像的重構(gòu)。通過(guò)計(jì)算光譜加權(quán)最小距離 (Spectral Angle Weighted Minimum Distance, SAWMD)，即按照公式(5～8)將h(q)與高光譜數(shù)據(jù)集的均值向量進(jìn)行匹配：

(5)

(6)

(7)

SAWMD(q)(i)=EMD(q)(i)·[1-cos(SAM)(q)(i)]n

(8)

式(5～8)中：EMD(q)(i)表示的是h(q)(i)中第i個(gè)重構(gòu)向量與均值向量的之間歐氏距離；SAM(q)(i)表示兩個(gè)向量的光譜角；參數(shù)n用于調(diào)節(jié)光譜角的權(quán)重，一般取值為1。光譜距離越小，表明結(jié)果越精確。按照公式(10)逐像元篩選轉(zhuǎn)換矩陣，實(shí)現(xiàn)高光譜影像的重構(gòu)，重構(gòu)高光譜擁有多光譜圖像的空間分辨率和幅寬。

(9)

SAWMD(q)(i)=min(SAWMD(1)(i),SAWMD(2)(i),…SAWMD(N)(i)

(10)

2.2 聚類后聚合的光譜分辨率增強(qiáng)方法:ISREM-SGA

在提取光譜數(shù)據(jù)集時(shí)，將高光譜圖像重采樣成多光譜影像大小，以提取對(duì)應(yīng)位置的像元光譜。但重采樣會(huì)導(dǎo)致像元發(fā)生變化，例如在使用最近鄰插值時(shí)，像元產(chǎn)生半個(gè)位置的偏移；使用雙線性內(nèi)插時(shí)也會(huì)破壞像元信息量，影響波譜分析[15]。在進(jìn)行空間配準(zhǔn)時(shí)，配準(zhǔn)影像的多源性與環(huán)境復(fù)雜性，也會(huì)導(dǎo)致較難獲得精確對(duì)準(zhǔn)的像元點(diǎn)[16]。為了改善上述不足，本研究提出一種基于光譜梯度角的改進(jìn)ISREM-SGA高光譜-多光譜圖像融合方法。具體操作如下：(1)按照SVM高光譜圖像分類結(jié)果，提取高光譜像元光譜，并根據(jù)光譜變換矩陣進(jìn)行降維。(2)計(jì)算降維后高光譜像元光譜，與多光譜影像上對(duì)應(yīng)位置及其鄰域內(nèi)像元光譜的梯度角來(lái)匹配同名點(diǎn)。計(jì)算光譜梯度角需對(duì)2個(gè)光譜向量先進(jìn)行一階求導(dǎo)，獲得其梯度向量，再計(jì)算2個(gè)梯度向量的廣義夾角[17]，光譜梯度角的計(jì)算如式(13)所示。該方法能夠反映光譜特征變化的優(yōu)勢(shì)，特別是對(duì)光譜曲線斜率的變化。將最小梯度角對(duì)應(yīng)的像元作為該高光譜像元最佳匹配點(diǎn)，提取出最具相似光譜特性的多光譜數(shù)據(jù)集。(3)按照上述公式(5～10)逐像元實(shí)現(xiàn)高光譜影像重構(gòu)。

SG(x)=(x2-x1,x3-x2,…，xn-x(n-1))

(11)

SG(y)=(y2-x1,y3-y2,…，yn-y(n-1))

(12)

(13)

3 結(jié)果與分析

3.1 影像視覺(jué)效果分析

融合圖像的質(zhì)量決定了圖像應(yīng)用的廣泛性，選取SREM和ISREM-SGA方法與常用的CNMF融合和Gram-Schmidt融合(以下簡(jiǎn)稱 GS融合)[18]進(jìn)行對(duì)比，其中CNMF與GS融合采用的是覆蓋范圍一致的高光譜和多光譜影像。彩色合成后的Chikusei(RGB分別對(duì)應(yīng)波段60(667 nm)、40(564 nm)、20(467 nm))和Pavia U(RGB分別對(duì)應(yīng)波段103(857 nm)、63(687 nm)、29(531 nm))融合圖像展示在下圖3中 。

圖 3 Chikusei與Pavia U彩色合成圖像Figure 3 Color composite images of Chikusei and Pavia U

從圖3中可以發(fā)現(xiàn)，ISRME-SGA融合影像的視覺(jué)效果與參考圖像比較接近，影像清晰，紋理特征與邊緣特征明顯，且色彩信息與參考圖像基本保持一致，不受空間變異度的影響。SREM方法生成Pavia U影像的植被和建筑物亮度值變低，發(fā)生了明顯的色彩畸變，而ISREM-SGA緩解了這一問(wèn)題。主要原因在于通過(guò)光譜梯度角提到的數(shù)據(jù)集，最具相似光譜特性，能夠有效地代表此類地物，從而獲得質(zhì)量更高的轉(zhuǎn)換矩陣用于影像重構(gòu)。ISREM-SGA方法與CNMF方法生成的融合影像，在視覺(jué)效果上與參考圖像最為相似，但I(xiàn)SREM-SGA在僅依靠重疊區(qū)域光譜和空間信息的基礎(chǔ)上，生成了高質(zhì)量融合影像，更具優(yōu)勢(shì)。在圖3中可以發(fā)現(xiàn)，ISRME-SGA融合影像在視覺(jué)效果上與參考圖像比較接近，影像清晰，紋理特征與邊緣特征明顯，且色彩信息與參考圖像基本保持一致，不受空間變異度的影響。SREM方法生成Pavia U影像植被和建筑物亮度值變低，發(fā)生了明顯的色彩畸變，而ISREM-SGA緩解了這一問(wèn)題。主要原因在于通過(guò)光譜梯度角提到的數(shù)據(jù)集，最具相似光譜特性，能夠有效地代表此類地物，從而獲得質(zhì)量更高的轉(zhuǎn)換矩陣用于影像重構(gòu)。ISREM-SGA方法與CNMF方法生成的融合影像，在視覺(jué)效果與參考圖像最為相似，但I(xiàn)SREM-SGA在僅依靠重疊區(qū)域光譜和空間信息的基礎(chǔ)上，生成了高質(zhì)量融合影像，更具優(yōu)勢(shì)。

3.2 典型地物光譜分析

分別提取4種融合影像的植被和建筑物光譜曲線，比較融合方法間的差異，光譜曲線如圖4所示。在圖4的光譜曲線對(duì)比中可以發(fā)現(xiàn)，在ISREM融合影像上提取的植被和建筑物光譜曲線，與參考光譜變化趨勢(shì)基本一致。從上圖2信噪比變化曲線中可以看出Chikusei影像在430～500 nm、700～730 nm，Pavia U影像在760～850 nm范圍內(nèi)高光譜波段的信噪比要小于其他波段，較低的光譜信息質(zhì)量導(dǎo)致光譜間產(chǎn)生差異。其中植被光譜相較于建筑物光譜誤差小，建筑物由于材質(zhì)不同，擁有不同的吸收或反射特性，對(duì)像元位置偏移敏感。在SREM生成的Pavia U融合影像中，受到像元偏移與復(fù)雜場(chǎng)景內(nèi)容的影響，純凈端元獲取困難，在進(jìn)行SVM分類時(shí)，影像無(wú)法識(shí)別建筑物純凈光譜，導(dǎo)致缺少相應(yīng)的光譜數(shù)據(jù)集，融合后建筑物光譜在除紅光之外的波段都發(fā)生明顯的畸變。盡管如此，ISREM方法的誤差值也要遠(yuǎn)小于SREM方法的誤差，由此說(shuō)明了ISREM方法對(duì)比原方法有了一定的改善。

圖 4 典型地物光譜曲線對(duì)比Figure 4 Comparison of spectrum of typical features

3.3 融合影像質(zhì)量統(tǒng)計(jì)分析

統(tǒng)計(jì)指標(biāo)上本研究選擇了相關(guān)系數(shù)R[19]、通用圖像質(zhì)量指標(biāo)(Universal Image Quality Index, UIQI)[20]、相對(duì)全局無(wú)量綱誤差(Erreur Relative Global Adimensionnelle Synthese，ERGAS)[21]、結(jié)構(gòu)相似度(Structural Similarity Index，SSIM)[22]、光譜角(Spectral Angle Mapper ,SAM)[23]對(duì)融合影像進(jìn)行質(zhì)量評(píng)價(jià)。一般認(rèn)為R、UIQI、SSIM值越接近1圖像失真越少，ERGAS和SAM則是值越小融合質(zhì)量越高。公式如下，計(jì)算結(jié)果見(jiàn)表2。

表 2 融合圖像質(zhì)量評(píng)價(jià)Table 2 Evaluation of fusion image

(14)

(15)

(16)

(17)

(18)

式(14～18)中：x表示融合圖像，y表示原始高光譜數(shù)據(jù)，x與y表示波段均值，σxy表示協(xié)方差，σx與σy表示方差，表示光譜向量的內(nèi)積，‖x‖、‖y‖表示光譜向量的范數(shù)，C1、C2是避免分母為0的常數(shù)，h、l表示多光譜圖像與高光譜圖像的空間分辨率，RMSEK表示均方根誤差，μK表示第K波段參考圖像的均值。

從表2中可以看出ISREM-SGA對(duì)比原方法，在空間和光譜指標(biāo)上都有所提升，其中Chikusei影像質(zhì)量略有提升，Pavia U影像在空間和光譜上的提升較明顯，R、UIQI、SSIM值分別增加了0.237、0.297、0.238，ERGAS、SAM值分別降低了35.81和7.312，具有較強(qiáng)的空間和光譜信息保持能力。ISREM-SGA對(duì)比CNMF融合差距較小，但I(xiàn)SREM-SGA方法只使用了部分重疊影像，運(yùn)算時(shí)間更短。相較于GS融合，ISREM-SGA在兩組數(shù)據(jù)集上都取得了較好的融合效果，更具穩(wěn)健性。選取融合前后每組對(duì)應(yīng)波段內(nèi)全部像元計(jì)算相關(guān)系數(shù)與信息熵(圖5)。從圖5中可以看出ISREM-SGA融合結(jié)果雖然在個(gè)別波段的值較低，但整體上有所提升，如450～700 nm、760～830 nm處，相關(guān)系數(shù)均在0.90以上。在相關(guān)系數(shù)差距大的波段，如900～950 nm處，信息熵也有明顯的差異，雖然熵值較大，但承載的是離散信息，在實(shí)際應(yīng)用中可以去除這些異常波段，以擴(kuò)大應(yīng)用范圍。

圖 5 相關(guān)系數(shù)與信息熵對(duì)比Figure 5 Correlation coefficient and information entropy

為進(jìn)一步比較各方法的光譜信息保持能力，進(jìn)行光譜角制圖[24]， 圖像中像元值代表光譜角大小，顏色越明亮光譜角越大(圖6)。通過(guò)光譜角制圖結(jié)果可以看出ISREM-SGA對(duì)比SREM方法，影像中框選的地物混合分布區(qū)和建筑區(qū)的光譜角降低，可較好地保持光譜信息。對(duì)比CNMF融合，ISREM-SGA的方法略有遜色，主要原因在于CNMF方法無(wú)須考慮重采樣和空間配準(zhǔn)造成的像元位置偏差，對(duì)像元分解的過(guò)程以及波譜分析的影響。

圖 6 各融合方法光譜角制圖 Figure 6 SAM map of different methods

為了衡量融合方法在空間信息上的保持能力，對(duì)融合結(jié)果進(jìn)行主成分變換。由于融合前后的影像在第1、2主分量中信息量都比較豐富，較難驗(yàn)證融合后增加的空間信息量，因此在圖7、8中展示變換后的第3、4主分量。根據(jù)主成分分量顯示圖，融合前影像的分量雖然可以看出部分紋理特征，但影像較模糊。SREM與ISREM-SGA融合方法雖然只使用重疊區(qū)域的空間信息，但對(duì)比其他清晰度增加，生成的Pavia U融合圖像，在圖8中框選的區(qū)域可以看出明顯的改善，表明了ISREM-SGA具有較好的空間信息保持能力，相較于其他融合方法有著一定的優(yōu)勢(shì)。

4 結(jié)論

采用ISREM-SGA融合、SREM融合、CNMF融合、GS融合進(jìn)行高光譜-多光譜的融合實(shí)驗(yàn)，并對(duì)融合圖像進(jìn)行光譜和空間信息質(zhì)量評(píng)價(jià)。評(píng)價(jià)結(jié)果表明：基于光譜梯度角的ISREM-SGA融合方法，光譜信息和空間信息保持能力較好，具有很好的穩(wěn)健性。與CNMF融合方法相比，融合圖像效果相似但運(yùn)算時(shí)間更短。ISREM-SGA方法通過(guò)光譜梯度角對(duì)像元光譜進(jìn)行匹配，能夠準(zhǔn)確地提取對(duì)應(yīng)位置的光譜數(shù)據(jù)集進(jìn)行融合研究，緩解了像元位置偏移造成的光譜畸變問(wèn)題，提高了融合圖像質(zhì)量。

但本方法仍存在不足，如重采樣的方法只選擇了最近鄰法，沒(méi)有對(duì)比影像在雙線性內(nèi)插和三次卷積內(nèi)插時(shí)融合結(jié)果的質(zhì)量。在今后的研究中，仍將進(jìn)一步改進(jìn)。

圖 7 Pavia U主成分變換后第3、4分量Figure 7 PC3,PC4 of Pavia U

圖 8 Chikusei主成分變換后第3、4分量Figure 8 PC3,PC4 of Chikuse