李昌元，劉國棟，譚 博

（1. 重慶交通大學(xué)土木工程學(xué)院，重慶 400074）

高光譜遙感影像在反映地物空間位置信息的同時(shí)也記錄著地物反射的光譜波段信息，兩種圖像信息相互結(jié)合，能更精細(xì)地表現(xiàn)光譜維中不同地物間的細(xì)微差異[1]，大幅提升地物的識(shí)別與分類精度。然而，高光譜遙感影像數(shù)據(jù)光譜波段眾多、特征空間維數(shù)過高，極易造成數(shù)據(jù)冗余，且光譜相鄰波段間存在較強(qiáng)的相關(guān)性，更會(huì)引發(fā)“維度災(zāi)難”[2]，給后續(xù)的地物識(shí)別與分類等帶來諸多困難。對(duì)高光譜遙感影像數(shù)據(jù)進(jìn)行降維，既能充分利用其豐富的數(shù)據(jù)量，又能降低數(shù)據(jù)處理代價(jià)、解決特征空間維數(shù)過高的問題，還能保留必要的地物信息，提高地物識(shí)別、分類的精度和效率。

高光譜遙感影像數(shù)據(jù)的降維技術(shù)主要包括波段選擇和波段提取兩種，按照不同的標(biāo)準(zhǔn)又可分為線性映射和非線性映射降維[2]?；诰€性變換降維的方法稱為線性映射降維，主要包括PCA[3]、線性判別分析（LDA）[4]和獨(dú)立成分分析（ICA）[5]等，方法簡(jiǎn)便且計(jì)算速度快。高光譜遙感數(shù)據(jù)本質(zhì)是非線性的，線性映射降維方法無法利用高光譜數(shù)據(jù)自身的非線性結(jié)構(gòu)。因此，諸多學(xué)者提出不同的非線性映射降維方法，主要包括基于核技術(shù)的非線性降維方法（如Fauvel M[6]等提出的KPCA、楊國鵬[7]等提出的核Fisher 判別分析）、DONG G J[8]等提出的等距映射（ISOMAP）、羅琴[9]等提出的局部線性嵌入（LLE）等。非線性映射降維方法彌補(bǔ)了線性映射降維方法不能發(fā)現(xiàn)高光譜影像數(shù)據(jù)非線性結(jié)構(gòu)的缺陷，但其處理復(fù)雜的非線性特征關(guān)系時(shí)計(jì)算量非常大，大幅增加了計(jì)算成本。

現(xiàn)有研究中，沒有明確對(duì)PCA 算法和KPCA 算法差異性的比較。本文通過設(shè)定閾值，計(jì)算了主成分分量的個(gè)數(shù)和降維時(shí)間，比較了兩種算法的數(shù)據(jù)壓縮能力、計(jì)算機(jī)運(yùn)行成本等內(nèi)部差異性；再加入多層次感知器（MLP）對(duì)降維后的數(shù)據(jù)進(jìn)行外部差異性分析，基于兩種降維算法的分類結(jié)果做進(jìn)一步評(píng)價(jià)，分析了兩種降維算法的優(yōu)劣。

1 基本原理與實(shí)現(xiàn)步驟

1.1 PCA算法

1.1.1 基本原理

PCA 算法是一種基于線性映射的特征提取技術(shù)。通過一定變換將高維影像數(shù)據(jù)變換到一個(gè)新的低維空間，使高維數(shù)據(jù)的最大方差投影在第一個(gè)低維空間的坐標(biāo)（即第一主成分分量）上，第二大方差投影在第二個(gè)低維空間的坐標(biāo)（第二主成分分量）上，以此類推[10]。利用少數(shù)幾個(gè)主成分分量將原始高維影像數(shù)據(jù)最大限度地保留下來，第一主成分分量包含了原始高維影像數(shù)據(jù)中的絕大部分信息[11]。PCA 算法主要利用協(xié)方差矩陣是一個(gè)實(shí)對(duì)角矩陣的性質(zhì)，即方差最大化、協(xié)方差最小化，來進(jìn)行降維。

1.1.2 實(shí)現(xiàn)步驟

假設(shè)高光譜遙感數(shù)據(jù)有P個(gè)波段，數(shù)據(jù)表示為X=(x1,x2,…,xN)=(X1,X2,…,XP)T，其中Xk為一個(gè)N×1 維的列向量（N為高光譜數(shù)據(jù)的像元數(shù)目，即N=m×n，m為行數(shù)，n為列數(shù)），則X為一個(gè)N×P的矩陣，降維后的低維輸出維度為d（d?P）。其主要步驟為[12-14]：

1）計(jì)算所有訓(xùn)練樣本的均值向量。

式中，i=1,2,…,N；k=1,2,…,P。

2）計(jì)算原始高維數(shù)據(jù)零均值化的標(biāo)準(zhǔn)矩陣。

3）計(jì)算所有訓(xùn)練樣本的協(xié)方差矩陣。

4） 計(jì)算協(xié)方差矩陣的特征值λ1,λ2,…,λP,λ1?λ2…?λP及其對(duì)應(yīng)的特征向量α1,α2,…,αP，并將特征值由大到小進(jìn)行排序，其特征向量會(huì)隨特征值的排序依次排列；通過得到的特征值，計(jì)算每個(gè)主成分所含的貢獻(xiàn)率σi和累計(jì)貢獻(xiàn)率δ。

5）取最大的d個(gè)特征值（d?P），將對(duì)應(yīng)的特征向量α1,α2,…,αd組成轉(zhuǎn)換矩陣A=[]α1,α2,…,αd，利用式（5）計(jì)算得到原始高光譜遙感影像數(shù)據(jù)X降至d維的數(shù)據(jù)Y。

變換后原始影像數(shù)據(jù)的絕大部分信息排在前面的幾個(gè)主成分分量中，眾多靠后分量包含的信息基本為噪聲[15]。因此，PCA 算法在一定程度上起到了降噪的作用。PCA算法的降維流程如圖1所示。

圖1 PCA算法降維流程圖

1.2 KPCA算法

1.2.1 基本原理

PCA算法一般適用于數(shù)據(jù)的線性降維，但高光譜數(shù)據(jù)的本質(zhì)是非線性的，線性映射降維方法無法表示其非線性結(jié)構(gòu)。為了實(shí)現(xiàn)高光譜數(shù)據(jù)的非線性特征提取，在線性映射降維方法的基礎(chǔ)上加入了核函數(shù)，將傳統(tǒng)線性映射降維方法核化。核化方法采用非線性映射將輸入空間的原始數(shù)據(jù)映射到高維特征空間，在特征空間中進(jìn)行線性處理，其關(guān)鍵在于引入核函數(shù)，忽略具體映射函數(shù)的表達(dá)式，直接得到低維數(shù)據(jù)映射到高維后的內(nèi)積（協(xié)方差矩陣的每個(gè)元素都是向量的內(nèi)積）[16]。目前，廣泛使用的核函數(shù)主要包括[17]：

1）線性核函數(shù)：K(x,y)=xTy+c。

2）多項(xiàng)式核函數(shù)：K(x,y)=[a(x·y)+c]d，c≥0,d≤N,a為常數(shù)。

4） 多層感知器核函數(shù)（Sigmoid 核函數(shù)）：K(x,y)=tanh[s(x·y)+c] ，s、c為常數(shù)。

KPCA算法是將輸入的矩陣數(shù)據(jù)X由一個(gè)非線性映射Φ 映射到更高維的特征空間? 中，使其線性可分，然后在高維特征空間? 上執(zhí)行線性PCA，利用核函數(shù)求得內(nèi)積，得到矩陣數(shù)據(jù)X在特征向量上的投影。其降維基本原理為[16-20]：

假設(shè)原始高光譜遙感影像數(shù)據(jù)矩陣為X=[x1,x2,…,xN]，其中每個(gè)樣本點(diǎn)Xi為P維列向量，X中共有N個(gè)樣本。利用一個(gè)非線性映射Φ 將矩陣X中的向量Xi映射到高維特征空間? （假設(shè)有K維），則有：

將矩陣X中所有樣本都映射到高維特征空間?中，得到一個(gè)K×N的新矩陣Φ(X)，即

假設(shè)矩陣Φ(X)已作中心化處理，則有：

在特征空間? 中，其協(xié)方差矩陣為：

式中，D?為K×K的矩陣。

利用式（8）求解特征空間? 的特征值問題：式中，K維列向量Ρ是特征空間?中對(duì)應(yīng)的特征向量。

將式（7）代入式（8），則有：

當(dāng)λ≠0 時(shí)，特征向量Ρ的線性表示為：

式中，α為N維列向量[α1,α2,…,αN]T。

此時(shí)，求解特征向量的問題隨之轉(zhuǎn)化為求解對(duì)應(yīng)α值的問題，將式（10）代入式（9），則有：

Φ(X)[Φ(X)]TΦ(X)α=λΦ(X)α

將兩邊左乘矩陣[Φ(X)]T，即

定義一個(gè)N×N的矩陣K=[Φ(X)]TΦ(X)，其第i行j列的元素為Ki,j=[Φ(xi)]TΦ(xj)，則式（12）可轉(zhuǎn)化為：

此時(shí)得到的特征向量α并不是PCA算法特征值分解時(shí)得到的基向量，而是利用特征空間映射樣本線性表示時(shí)的線性系數(shù)。實(shí)際上應(yīng)求得α對(duì)應(yīng)的特征向量P的標(biāo)準(zhǔn)向量，因此需對(duì)Ρ進(jìn)行歸一化處理，即

PTP=1

則有：

求解得到K的特征值和特征向量α后，再利用式（10）求解得到特征向量Ρ，在KPCA算法特征空間? 中的主成分方向就是Ρ，則原始數(shù)據(jù)映射到特征空間? 上第j維向量Ρj的投影為：

若特征空間? 中的數(shù)據(jù)未滿足中心化條件，則需對(duì)矩陣K進(jìn)行改正，將K?代替式（13）中的K，即

可轉(zhuǎn)化為：

式中，I1N為一個(gè)N×N的矩陣（與核矩陣維數(shù)相同），其所有元素均為1N。

1.2.2 實(shí)現(xiàn)步驟

若高光譜遙感影像數(shù)據(jù)有M個(gè)訓(xùn)練樣本，每個(gè)樣本有P個(gè)波段，記X=[X1,X2,…,XM]，將其表示為M×P的矩陣，降維后的低維輸出維度為d。其主要步驟為[21]：

1） 將矩陣X映射到特征空間? 中，即Φ(X)=[Φ(x1),Φ(x2),…,Φ(xM)]。

2）選擇合適的核函數(shù)，對(duì)核矩陣作中心化處理，建立標(biāo)準(zhǔn)核矩陣，并計(jì)算標(biāo)準(zhǔn)核矩陣，即=K-I1M K-KI1M+I1M KI1M。

3）計(jì)算核矩陣K?的特征值λ1,λ2,…,λP,λ1λ2…?λP及其對(duì)應(yīng)的特征向量Ρ1,Ρ2,…,ΡP，并將得到的特征值由大到小進(jìn)行排序，其特征向量會(huì)隨特征值的排序依次排列。

4） 對(duì)特征向量Ρ1,Ρ2,…,ΡP進(jìn)行正交單位化（Gram-Schmidt 法），得到的特征向量即為核主成分分量α1,α2,…,αP。

5）取最大的d個(gè)特征值（d?P），提取對(duì)應(yīng)的d個(gè)主成分分量α1,α2,…,αd，求得投影矩陣Y，即Y=，α=(α1,α2,…,αd)，Y為矩陣X數(shù)據(jù)經(jīng)KP?CA算法降維后的數(shù)據(jù)。

KPCA算法降維流程如圖2所示。

圖2 KPCA算法降維流程圖

2 實(shí)驗(yàn)方法與結(jié)果分析

2.1 實(shí)驗(yàn)數(shù)據(jù)

本文選用1992 年AVIRIS 拍攝的Indian Pines（印度松樹）數(shù)據(jù)集。該影像數(shù)據(jù)集截取145像素×145像素大小作為高光譜遙感影像數(shù)據(jù)降維的測(cè)試數(shù)據(jù)，共有21 025個(gè)像素，其中背景像素10 776個(gè)，地物像素10 249個(gè)；共有220個(gè)波段，但其中20個(gè)波段受噪聲和水汽影響，因此剔除這20個(gè)波段，選用剩余的200個(gè)波段作為研究對(duì)象。Band43、Band23 和Band10 三個(gè)波段疊加形成的假彩色合成圖像如圖3a所示；印度松樹數(shù)據(jù)集16種地物的真實(shí)地物圖像如圖3b所示。

圖3 印度松樹測(cè)試數(shù)據(jù)集

2.2 實(shí)驗(yàn)結(jié)果分析

實(shí)驗(yàn)利用Python 編寫PCA 和KPCA 兩種降維算法，利用印度松樹數(shù)據(jù)集進(jìn)行測(cè)試。為了比較PCA算法與KPCA 算法之間的差異性，設(shè)置固定的閾值T（即累積貢獻(xiàn)率）≥0.95，得到降維后各主成分分量的貢獻(xiàn)率和累計(jì)貢獻(xiàn)率，重復(fù)5 次得其平均值，結(jié)果如表1所示。

表1 PCA算法和KPCA算法降維后各主成分分量的貢獻(xiàn)率與累計(jì)貢獻(xiàn)率

由表1 可知，原印度松樹影像數(shù)據(jù)200 維包含的信息，經(jīng)PCA算法降維后其數(shù)據(jù)前5個(gè)主成分分量就表達(dá)了整個(gè)影像數(shù)據(jù)95.038 1%的信息，而經(jīng)KPCA算法降維后其數(shù)據(jù)僅前3 個(gè)主成分分量就包含原始影像數(shù)據(jù)95.286 9%的信息量，說明PCA算法和KPCA算法均能對(duì)高光譜遙感影像數(shù)據(jù)進(jìn)行優(yōu)化處理，達(dá)到數(shù)據(jù)壓縮、約簡(jiǎn)的目的，而KPCA 算法作為PCA 算法的核函數(shù)延伸方法，處理高維影像數(shù)據(jù)時(shí)在盡可能包含更多原始信息的基礎(chǔ)上具有更強(qiáng)的壓縮能力和的更好降維效果。

通過計(jì)算兩種降維算法的時(shí)間，比較了二者的時(shí)間復(fù)雜度。兩種算法均運(yùn)行5次取平均值，對(duì)比其運(yùn)行成本，PCA算法用時(shí)16 s，KPCA算法用時(shí)85.6 s，PCA算法的優(yōu)勢(shì)較明顯。原因在于，KPCA算法做了非線性變換，在高維空間進(jìn)行協(xié)方差矩陣的特征值分解，再采用與PCA 算法相同的方式進(jìn)行降維，在高維空間中KP?CA算法運(yùn)用了核函數(shù)運(yùn)算，計(jì)算量遠(yuǎn)超于PCA算法。

為進(jìn)一步對(duì)比兩種降維算法的差異性，充分分析二者的特征差異性，本文選擇基于分類應(yīng)用的降維效果評(píng)估方法，以MLP為分類器進(jìn)行分類。MLP分類器是目前遙感領(lǐng)域的主流分類器之一，由許多相同的處理單元并聯(lián)組合而成，能進(jìn)行大量簡(jiǎn)單單元的并行活動(dòng)，對(duì)信息的處理能力、普適性較強(qiáng)，較適合高光譜遙感數(shù)據(jù)這類大型數(shù)據(jù)的分類處理，在遙感影像分類中有較多的應(yīng)用實(shí)例，更容易驗(yàn)證PCA 算法和KPCA算法的降維效果。

兩種降維算法基于MLP分類應(yīng)用的降維效果如圖4所示，可以看出，第1、7和8種地物的分類準(zhǔn)確率為100%，由于這3 種地物樣本太少，無法做到有效預(yù)測(cè)，因此其分類沒有價(jià)值。地物序號(hào)為0 的是影像的背景值，是分類中錯(cuò)分的主要因素。總體而言，KP?CA算法經(jīng)分類后各地物的分類精度略優(yōu)于PCA算法。

圖4 基于MLP分類應(yīng)用的降維效果對(duì)比

由兩種降維算法分別經(jīng)MLP分類器分類后的結(jié)果可知，KPCA 算法的總體精度高于PCA 算法，錯(cuò)誤分類比例也較少，但分類的時(shí)間略長(zhǎng)，如表2所示。

表2 基于兩種降維數(shù)據(jù)的MLP分類結(jié)果精度統(tǒng)計(jì)表

3 結(jié) 語

1）PCA算法和KPCA算法通過舍去一部分原始影像中不必要的數(shù)據(jù)信息來降低維度，使得樣本采樣更密集，緩解了維度災(zāi)難，降低了噪聲影響。實(shí)驗(yàn)結(jié)果表明，在相同數(shù)量的主成分分量條件下，KPCA 算法能更很好地包含高光譜遙感影像中的原始信息，具有比PCA算法更好的信息集中性。

2）KPCA算法在低維空間應(yīng)用核函數(shù)，巧妙地從低維升至高維，在高維空間中實(shí)現(xiàn)非線性映射，但大幅增加了算法的復(fù)雜度，降低了高光譜遙感影像數(shù)據(jù)降維處理的時(shí)效性。在這一點(diǎn)上，PCA算法計(jì)算量小且速度快，具有更大的優(yōu)勢(shì)。

3）KPCA算法提取的各分量特征在后續(xù)的影像分類中具有良好的表現(xiàn)，各地物總體精度均高于PCA算法，具有足夠的適用性。

4）PCA 算法在提取遙感數(shù)據(jù)主要信息的同時(shí)將所有樣本當(dāng)作一個(gè)整體看待，尋找均方誤差最小的線性映射投影，卻忽略了地物的類別屬性，導(dǎo)致后續(xù)的分類精度降低，出現(xiàn)比KPCA算法更多的錯(cuò)分和漏分現(xiàn)象。

由于高光譜數(shù)據(jù)特征維度過高，對(duì)后續(xù)的研究易造成巨大影響，因此找到一個(gè)既能充分利用高光譜豐富數(shù)據(jù)量又能降低消耗成本的降維方法才能解決處理高光譜數(shù)據(jù)成本大的問題。同時(shí)，選擇不同的核函數(shù)也會(huì)不同程度地影響降維與分類的效果。因此，非線性KPCA算法的參數(shù)選擇也是一個(gè)降維算法的改進(jìn)方向。