基于PCA 和KPCA 的高光譜遙感數據降維對比研究

李昌元，劉國棟，譚 博

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

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

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

現有研究中，沒有明確對PCA 算法和KPCA 算法差異性的比較。本文通過設定閾值，計算了主成分分量的個數和降維時間，比較了兩種算法的數據壓縮能力、計算機運行成本等內部差異性；再加入多層次感知器（MLP）對降維后的數據進行外部差異性分析，基于兩種降維算法的分類結果做進一步評價，分析了兩種降維算法的優劣。

1 基本原理與實現步驟

1.1 PCA算法

1.1.1 基本原理

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

1.1.2 實現步驟

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

1）計算所有訓練樣本的均值向量。

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

2）計算原始高維數據零均值化的標準矩陣。

3）計算所有訓練樣本的協方差矩陣。

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

5）取最大的d個特征值（d?P），將對應的特征向量α1,α2,…,αd組成轉換矩陣A=[]α1,α2,…,αd，利用式（5）計算得到原始高光譜遙感影像數據X降至d維的數據Y。

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

圖1 PCA算法降維流程圖

1.2 KPCA算法

1.2.1 基本原理

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

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

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

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

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

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

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

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

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

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

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

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

當λ≠0 時，特征向量Ρ的線性表示為：

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

此時，求解特征向量的問題隨之轉化為求解對應α值的問題，將式（10）代入式（9），則有：

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

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

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

此時得到的特征向量α并不是PCA算法特征值分解時得到的基向量，而是利用特征空間映射樣本線性表示時的線性系數。實際上應求得α對應的特征向量P的標準向量，因此需對Ρ進行歸一化處理，即

PTP=1

則有：

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

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

可轉化為：

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

1.2.2 實現步驟

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

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

2）選擇合適的核函數，對核矩陣作中心化處理，建立標準核矩陣，并計算標準核矩陣，即=K-I1M K-KI1M+I1M KI1M。

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

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

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

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

圖2 KPCA算法降維流程圖

2 實驗方法與結果分析

2.1 實驗數據

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

圖3 印度松樹測試數據集

2.2 實驗結果分析

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

表1 PCA算法和KPCA算法降維后各主成分分量的貢獻率與累計貢獻率

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

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

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

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

圖4 基于MLP分類應用的降維效果對比

由兩種降維算法分別經MLP分類器分類后的結果可知，KPCA 算法的總體精度高于PCA 算法，錯誤分類比例也較少，但分類的時間略長，如表2所示。

表2 基于兩種降維數據的MLP分類結果精度統計表

3 結 語

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

2）KPCA算法在低維空間應用核函數，巧妙地從低維升至高維，在高維空間中實現非線性映射，但大幅增加了算法的復雜度，降低了高光譜遙感影像數據降維處理的時效性。在這一點上，PCA算法計算量小且速度快，具有更大的優勢。

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

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

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