相對熵結(jié)合互信息的高光譜波段選擇方法

袁博，胡冰

(1.南陽理工學(xué)院 計算機(jī)與信息工程學(xué)院，河南 南陽 473004；2.南陽理工學(xué)院 學(xué)報編輯部，河南 南陽 473004)

0 引言

研究發(fā)現(xiàn)，光譜往往在很大程度上表征了地物的本征特性[1]，光譜分辨率的提高有助于對地物的精確識別和分類[2]。然而，高光譜遙感數(shù)據(jù)的數(shù)據(jù)量大、冗余信息多，高光譜影像處理算法普遍具有約束條件多、運(yùn)算量大、計算復(fù)雜等特點，為了降低計算復(fù)雜度并改善處理精度，不僅需要從原始數(shù)據(jù)中提取地物光譜信息，還應(yīng)該在保留原始光譜結(jié)構(gòu)與有效光譜信息的前提下，盡可能減少參與運(yùn)算的數(shù)據(jù)量，降低光譜信息冗余[3]。

高光譜的數(shù)據(jù)降維方法主要可分為維數(shù)約簡和波段選擇兩大類。維數(shù)約簡就是通過線性或非線性數(shù)學(xué)變換，將高維數(shù)據(jù)轉(zhuǎn)化為可解的低維數(shù)據(jù)[4]。波段選擇則是從多達(dá)數(shù)百個波段的高光譜影像中，選擇具有較好分類識別能力的波段組合，代替原始全波段數(shù)據(jù)參與分析與應(yīng)用[5]。

維數(shù)約簡通常涉及數(shù)學(xué)運(yùn)算，容易破壞數(shù)據(jù)原始光譜結(jié)構(gòu)，不利于地物光譜信息提取，波段選擇則不存在此問題。波段選擇方法的研究成果很多，包括基于稀疏的波段選擇[6](sparse based band selection，SpaBS)、稀疏矩陣分解的波段選擇法[7](sparse nonnegative matrix factorization，SNMF)、自適應(yīng)波段選擇[8](adaptive band selection，ABS)、加權(quán)概率原型分析的高光譜影像波段選擇[9]等。

不同于上述方法，本文提出一種相對熵結(jié)合互信息的光譜綜合信息量分析波段選擇(Kullback-Leibler divergence and mutual information based band selection，KLMIBS)方法來研究高光譜影像的波段選擇問題。通過構(gòu)建高光譜影像光譜信息的總體分布模型進(jìn)行波段選擇，并利用一個高光譜影像公開數(shù)據(jù)集進(jìn)行實驗與精度分析，驗證本文提出的KLMIBS方法的有效性。

1 光譜信息分布模型

高光譜影像波段選擇的結(jié)果中，除了要求不改變原始光譜數(shù)據(jù)結(jié)構(gòu)，保持原始波段物理含義不變之外，還應(yīng)該盡可能多地保留地物的有效光譜信息[10]。因此，需要對整幅影像的光譜信息總體分布情況進(jìn)行分析與建模。KLMIBS方法選擇相對熵(relative entropy，RE)[11]結(jié)合互信息[12](mutual information，MI)構(gòu)建光譜信息總體分布的定量描述模型，指導(dǎo)高光譜影像的波段選擇。

1.1 待選擇的波段數(shù)量

研究波段選擇問題，首先需要確定待選擇的波段數(shù)量。高光譜影像處理算法通常要求選出的波段數(shù)量必須大于或等于端元數(shù)量[13]。端元數(shù)量可通過主成分分析[14](principal component analysis，PCA)、最小噪聲分離[15](minimum noise fraction，MNF)、奇異值分解[16](singular value decomposition，SVD)、噪聲子空間投影[17](noise subspace projection，NSP)和基于最小誤差高光譜信號識別[18](hyperspectral signal identification by minimum error，HySime)等方法得到。

其中，HFC和NWHFC方法都是基于相關(guān)矩陣和協(xié)方差矩陣特征值之間的協(xié)方差趨近于0的假設(shè)條件，利用虛擬維度估計端元數(shù)量的。不同之處在于，NWHFC在執(zhí)行HFC之前，對樣本相關(guān)矩陣和樣本協(xié)方差矩陣進(jìn)行了噪聲白化，去除了噪聲方差的相關(guān)性，使得虛擬維度的估計更加精確。HFC和NWHFC的主要問題在于，估計的準(zhǔn)確度比較一般，且樣本數(shù)量不是很大時，其假設(shè)條件不成立。NSP僅采用白化后的樣本協(xié)方差矩陣進(jìn)行虛擬維數(shù)的估計，不受樣本數(shù)量的影響。在噪聲估計正確的前提下，尤其是樣本個數(shù)不多時，該方法優(yōu)于HFC和NWHFC。

對于高光譜遙感影像的端元數(shù)量估計問題，PCA丟棄的主分量中可能包含真實端元；HFC和NWHFC的假設(shè)條件受樣本數(shù)量影響很大；NSP 則需要精確估計噪聲和給定誤警率。HySime盡管計算過程較復(fù)雜，但不需要任何參數(shù)且估計準(zhǔn)確度較高，綜合性能最優(yōu)。因此，選擇HySime算法估計高光譜影像的端元數(shù)量。

1.2 高光譜影像光譜差異信息模型

相對熵又稱KL散度(Kullback-Leibler divergence)或信息散度(information divergence)，是2個概率分布P和Q差別的非對稱性的度量，用來度量使用基于Q的編碼來編碼來自P的樣本，平均所需的額外比特個數(shù)。為避免混淆，下文將相對熵、KL散度和信息散度統(tǒng)稱為KL散度。

設(shè)P(x)和Q(x)是X取值的2個離散概率分布，則P對Q的KL散度為：

(1)

KL散度僅當(dāng)概率P和Q各自總和均為1，且對于任何x都滿足Q(x)>0及P(x)>0時，才有定義。

對于連續(xù)隨機(jī)變量，其概率分布P和Q可按積分方式定義為式(2)：

(2)

式中：p和q分別表示分布P和Q的密度。

從式(1)和式(2)可以看出，KL散度總是非負(fù)，當(dāng)且僅當(dāng)分布P和Q完全相同時，KL散度為0；P和Q的相似度越高，KL散度越??；相似度越低，KL散度越大?？梢詫⒎植糚和Q分別與高光譜圖像中2個波段的光譜數(shù)據(jù)進(jìn)行類比，則2波段數(shù)據(jù)的相似度越低，信息差異(即有效信息量)越大，則對應(yīng)的KL散度數(shù)值越大，反之則越小。因此，選擇利用KL散度參數(shù)描述波段間的有效光譜信息。

DKL(P‖Q)只能描述2個波段P和Q之間的光譜信息差異，而波段選擇過程需要對高光譜影像全部波段間光譜信息差異的整體分布進(jìn)行定量描述。因此，定義了一個KL散度矩陣MKL。對于高光譜數(shù)據(jù)X∈RN×M(R代表非負(fù)實數(shù)集，N是全部像元的數(shù)量，M是所有波段的數(shù)量)，MKL是一個M×M的方陣，如式(3)所示。

(3)

式中:Xi表示高光譜數(shù)據(jù)中所有像元在第i波段的數(shù)據(jù)；DKL(Xi‖Xj)表示所有像元從波段i到波段j的KL散度，可用公式(1)計算。

1.3 高光譜影像光譜冗余信息模型

為了更全面地描述光譜信息量整體分布情況，不僅需要建立有效光譜信息的整體分布模型，還需要對冗余光譜信息進(jìn)行整體分布的定量化描述，建立光譜冗余信息的整體分布定量模型。為此，論文引入了信息學(xué)中的互信息MI(mutual information)參數(shù)。

互信息是用來度量隨機(jī)變量之間獨立性的基本準(zhǔn)則，可以表示成KL散度的形式。多個隨機(jī)變量之間的互信息定義為其聯(lián)合概率密度函數(shù)與各邊緣密度函數(shù)乘積之間的KL散度，即：

(4)

將式(4)右邊繼續(xù)展開：

(5)

由于KL散度的非負(fù)性可知：

(6)

互信息量I(ai；bj)在聯(lián)合概率空間P(XY)中的統(tǒng)計平均值，稱為平均互信息(average mutual information，AMI)，其計算公式如(7)所示。

(7)

利用平均互信息量是為了克服互信息量的隨機(jī)性，使其成為一個確定的量。關(guān)于平均互信息量的物理意義，以I(X；Y)為例，即在對Y一無所知的情況下，X的先驗不定度與收到Y(jié)后關(guān)于X的后驗不確定度之差?？梢杂闷骄バ畔⒘縼矶勘硎咀V間冗余信息，為波段選擇提供依據(jù)。

基于上述分析，本文定義了一個平均互信息矩陣，記為MI，作為高光譜影像整體光譜冗余信息模型。MI是一個M行、M列的方陣，如式(8)所示。

(8)

式中:I(Xi；Xj)表示波段i和波段j之間所有像元反射率的平均互信息，該數(shù)值可依據(jù)式(7)進(jìn)行計算。

1.4 高光譜影像光譜綜合信息模型

為了更加全面地表述所有波段之間的光譜信息分布，同時抑制單參數(shù)的負(fù)面影響，算法將KL散度矩陣MKL和平均互信息量矩陣MI以不同的權(quán)重系數(shù)進(jìn)行結(jié)合，生成一個綜合光譜信息矩陣S，作為高光譜影像的整體光譜信息模型，并將其作為波段選擇過程的主要參考。

通過對KL散度和平均互信息參數(shù)的綜合，可以完整地反映光譜信息量的譜間相關(guān)特征。式(9)給出了圖像中任意2個譜段數(shù)據(jù)X和Y之間的譜間相關(guān)模型。

S(X，Y)=DKL(P‖Q)-λ·MI(X；Y)

(9)

式中：S(X，Y)為波段X和波段Y之間的譜間相關(guān)性；DKL(P‖Q)為波段X的概率分布P和波段Y的概率分布Q之間的KL散度；MI(X；Y)為波段X和波段Y之間的平均互信息；λ為權(quán)重系數(shù)，由DKL(P‖Q)和MI(X；Y)的均值的相對幅度大小決定。

容易知道，與MKL和MI相同，S也是一個M行、M列的方陣，如式(10)所示。

(10)

S中較大的數(shù)值代表較大的光譜信息差異(也可以理解為有效光譜信息，區(qū)別于冗余光譜信息)和對應(yīng)波段之間較小的光譜信息依賴程度(也可以理解為較小的光譜信息冗余)。S(i，j)表示波段i和波段j之間所有像元的綜合光譜信息差異(或可以理解為有效光譜信息)。S(i，j)的計算如式(11)所示。

S(i，j)=cKLDKL(Xi‖Xj)-cIMI(Xi；Xj)

(11)

式中：cKL和cI分別對應(yīng)KL散度矩陣MKL和平均互信息量矩陣MI的權(quán)重系數(shù)。由于KL散度和平均互信息量對于波段選擇過程的參考意義基本上同等重要，為了使對應(yīng)矩陣MKL和MI對于S的貢獻(xiàn)度基本相當(dāng)，需要參照歸一化思想，對MKL和MI按照其內(nèi)部元素的平均幅值的比例關(guān)系分別賦予相應(yīng)權(quán)重。具體來說，即權(quán)重系數(shù)cKL和cI由MKL和MI中所有元素絕對值的均值的比值確定，如式(12)所示。

(12)

由于只關(guān)心cKL和cI之間的相對比值，而非二者的絕對數(shù)值，因此S的計算公式可以轉(zhuǎn)換為式(13)。

(13)

1.5 KLMIBS處理流程

希望利用矩陣S同時選出全部待選波段是不切實際的。本文算法主要基于如下技術(shù)思路：每次僅從剩余波段中選擇一個波段，該波段在所有待選波段中，包含了相對于已選擇波段集合最多的有效光譜信息。

對于矩陣S中的隨機(jī)一列：si={S(1，i)，S(2，i)，…，S(M，i)}，它代表波段i相對于其余全部波段的有效光譜信息(差異化光譜信息)。

S是MKL和MI的綜合，S中較大數(shù)值的元素，意味著對應(yīng)2個波段之間同時具有較大的KL散度值和較小的平均互信息，也就是說，對應(yīng)2個波段之間的光譜信息差異較大，光譜信息冗余較小，這正是波段選擇算法的理想依據(jù)。

KLMIBS波段選擇方法流程如下：

①利用Hysime估計高光譜圖像的端元數(shù)量，將端元數(shù)量估計結(jié)果確定為波段選擇數(shù)量T。

②計算綜合信息量矩陣S?；谑?1)和式(3)計算KL散度矩陣MKL，基于式(7)和式(8)計算平均互信息矩陣MI，基于式(12)和式(13)計算S。

③選擇第一個波段。選擇x1∈[1，M]，M為波段數(shù)量，使得sx1={S(1，x1)，…，S(M，x1)}在所有的si={S(1，i)，…，S(M，i)}，i=1，2，…，M中取最大值，則序號為x1的波段即為第一個被選擇波段。

④依次選擇剩余T-1個波段。假設(shè)全部已選擇波段的序號集合為K={K1，…，KU}，其中U為已選波段數(shù)量，選擇xi∈[1，M]，i∈[2，T]，xi?K，計算SK=S(K1，xi)+…+S(KU，xi)，使得SK在所有的xi中取最大值，則xi即為第i個被選擇波段序號。按照④中步驟依次選擇剩余的T-1個波段，并將選中的每一個波段序號加入集合K，則最終生成的波段序號集合K即為算法的波段選擇結(jié)果。

2 真實數(shù)據(jù)波段選擇實驗

真實實驗數(shù)據(jù)為高光譜數(shù)據(jù)圖像收集實驗儀器(hyperspectral digital imagery collection experiment，HYDICE)在美國華盛頓特區(qū)采集到的高光譜數(shù)據(jù)集，如圖1所示。該高光譜影像的行數(shù)與列數(shù)均為400。該數(shù)據(jù)原始波段數(shù)量210，波段范圍0.4～2.4 μm，波段寬度為10 nm，涵蓋了可見光和近紅外譜段范圍；在去除了0.9～1.4 μm范圍內(nèi)的大氣吸收波段后，剩余191個有效波段。

圖1 美國華盛頓特區(qū)HYDICE高光譜數(shù)據(jù)

根據(jù)KLMIBS波段選擇方法流程，首先利用HySime方法得到該影像子區(qū)的端元數(shù)量估計值為10，確定待選擇的波段數(shù)量為10。在191個有效波段中，KLMIBS最終選擇的10個波段的下標(biāo)集合如下：{5，13，26，39，47，56，71，82，136，191}。

圖2給出了波段選擇算法結(jié)果對應(yīng)的綜合光譜信息矩陣S的三維視圖(圖2(a))，以及10個選擇波段的光譜三維視圖(圖2(b))。如圖2(a)所示，圖中峰谷起伏程度越大，代表各波段之間的光譜差異信息(有效信息)越多，光譜相關(guān)信息(冗余信息)越少；說明KLMIBS方法中的矩陣S能夠?qū)蓛刹ǘ沃g的光譜信息差異有效凸顯出來，同時能夠?qū)λ胁ǘ沃g的光譜有效信息整體分布情況進(jìn)行定量描述。如圖2(b)所示，10個選定波段的集合中各波段之間的光譜幅值差異十分顯著，說明選出的波段集合包含的光譜信息冗余小，有效光譜信息量大。

圖2 KLMIBS波段選擇實驗結(jié)果

為進(jìn)一步定量分析和驗證KLMIBS方法的性能，將KLMIBS與前文所述的SpaBS、SNMF和ABS 3種代表性波段選擇方法的波段選擇結(jié)果在真實高光譜遙感圖像中的像元解混精度進(jìn)行對比。

像元解混算法統(tǒng)一采用最小體積約束的非負(fù)矩陣分解(minimum volume constrained NMF，MVCNMF)。MVCNMF是目前光譜解混研究中成果最多的基于凸面幾何學(xué)理論的方法中很有代表性的一種，通過將體積限制加入到NMF中，將最小二乘分析和凸面幾何結(jié)合起來，具有較高的解混精度。MVCNMF提出的代價函數(shù)包括兩部分，一部分為估量觀測數(shù)據(jù)與端元和豐度重建數(shù)據(jù)之間的近似誤差，另一部分由最小體積限制組成。

利用MVCNMF的端元光譜估計值與地面實測真值之間的光譜夾角距離[13](spectral angle distance，SAD)作為解混結(jié)果的精度評價參數(shù)，單位為弧度。SAD越小，代表像元解混精度越高，解混效果越好。

KLMIBS與SpaBS、SNMF和ABS的波段選擇結(jié)果的MVCNMF解混精度對比如表1所示。

表1 不同波段選擇算法的MVCNMF解混精度

如表1所示，與其他3種方法相比，KLMIBS對應(yīng)的SAD最小，代表其選擇結(jié)果的MVCNMF解混精度最優(yōu)，說明該方法能夠很好地保留有效光譜信息，去除冗余光譜信息，波段選擇效果優(yōu)于多數(shù)同類算法。

3 結(jié)束語

論文研究結(jié)果表明，在波段選擇過程中對高光譜影像進(jìn)行光譜信息量分析，通過引入KL散度與(平均)互信息2種信息學(xué)參數(shù)，分別對波段間的光譜差異信息與光譜冗余信息進(jìn)行建模，并生成反映整幅影像有效光譜信息整體分布情況的綜合模型，能夠在數(shù)量非常有限的波段集合中，相對于其余幾種代表性波段選擇方法保留更多的光譜有效信息，有助于后續(xù)高光譜影像分析與處理。

相對于前人在高光譜影像波段選擇方面的研究成果，論文主要在以下兩個方面進(jìn)行了改進(jìn)與創(chuàng)新：一是通過賦予不同權(quán)重系數(shù)，融合2種代表性信息學(xué)參數(shù)，對高光譜影像的有效光譜信息整體分布狀況進(jìn)行建模和定量描述；二是提出的KLMIBS方法屬于非監(jiān)督方法，能夠顯著降低對先驗知識的依賴程度，相對于部分對先驗知識依賴程度較高的波段選擇方法，在一定程度上改善了實用性，在后續(xù)處理中，其波段選擇結(jié)果也具有令人滿意的精度。

在今后工作中，該方法的適用范圍與穩(wěn)定性還需進(jìn)一步分析與驗證；基于KLMIBS的技術(shù)思路，其他信息學(xué)或統(tǒng)計學(xué)參數(shù)組合，以及其他權(quán)重系數(shù)計算方法，也都值得進(jìn)一步嘗試與檢驗。