基于空間統計學的高光譜遙感影像主成分選擇方法

孫肖， 彭軍還， 趙鋒， 王曉陽， 呂潔， 張登峰

(1.中國地質調查局廊坊自然資源綜合調查中心，廊坊 065000； 2.中國地質大學(北京)土地科學技術學院，北京 100083； 3.中國地質調查局烏魯木齊自然資源綜合調查中心，烏魯木齊 830057； 4.中國地質調查局西安礦產資源調查中心,西安 710100)

0 引言

在面向任務的高光譜遙感影像數據分析中，由于高光譜遙感影像波段數相對比較多，龐大的數據量給后續的分析處理帶來了極大的挑戰。在高光譜分類中，把分類精度隨著波段數量的增加先升高后降低的現象叫作“Hughes”現象[1]。Chang[2]研究發現最高有94%的波段是可以舍棄的，而且不會影響分類的精度。因此，在高光譜遙感影像的研究中，一般會首先進行降維處理，主成分分析(principal component analysis，PCA)就是一種常用的線性降維方法[3]。

PCA將數據的方差作為線性變換的標準，因此，一般按照變換后數據的累積方差進行主成分的選擇。降維主要的目的就是降低數據維數的同時，盡可能保留信息，而選擇方差較大的主成分必然會帶來信息的損失。目前沒有一種有效的方法來決定該選擇哪個主成分。

Jolliffe[4]通過大量實驗研究將PCA變換后主成分選擇的經驗閾值定為累計方差貢獻率大于0.85，但是該閾值在高光譜遙感的研究中具有局限性。PCA在高光譜遙感領域的應用主要包含兩大類： 一類是應用于解混、變化檢測、數據壓縮、目標探測、去噪等的研究中，根據研究的目的和內容不同，一般選擇某一特定主成分或某幾個主成分進行研究[5-10]； 另一類是應用于分類研究，Chang等[11]認為利用累計方差貢獻率大于0.99的主成分進行分類研究效果比較好。Li等[12]認為累計方差貢獻率大于0.9就能保證分類精度大于0.85。臧卓等[13]對高光譜遙感影像降維后的主成分進行分類測試，發現累計方差貢獻率與分類精度沒有必然聯系，而主成分的個數對分類結果的影響較為明顯，認為保留前15～20個主成分較為合適。黃鴻等[14-16]認為可以給定一定數量的主成分進行分類。臧卓等[17-19]逐次增加主成分數量進行分類，根據分類精度確定合適的主成分個數，該方法雖然能保證分類精度最高，但是效率較慢。Mather等[20]指出不能僅依靠特征值對應的主成分來做圖像分類，還應考慮圖像的實際視覺效果。Rodarmel等[21]分別采用編號1—5、1—10、1—25、1—50的主成分分段計算了分類的精度，認為可以用5%～10%的主成分個數進行分類。Ibarrola-Ulzurrun等[22]將常用的主成分選擇方法分為4類(基于特征值、紋理特征、類別變換和感興趣區分離)，分別利用前2，5，10，15，20主成分對這4類方法的分類精度進行了對比研究，認為特征值不是最適合的主成分選擇方法，類別變換和感興趣區分離需要人為的確定感興趣區，因此紋理特征是比較適合的主成分選擇方法，典型的紋理特征指標即信息熵。

以上主成分選擇方法雖然取得了一定的效果，但是仍存在依據不充分、效率較低、主觀性較強的問題。而且，PCA變換結果不隨噪聲排列，方差也不能判斷噪聲的大小。以上方法會導致部分圖像質量較好的主成分被舍去，而部分圖像質量較差的主成分參與分類的現象。從實際應用來看，編號較大的一些主成分對分類結果也有一定的影響[23]。

目前，定量的PCA變換后主成分選擇方法的研究還比較少。Zheng等[24]提出 GA-Fisher算法，能有效地增加編號較大的有效主成分。Zhang等[25]利用免疫克隆選擇算法對主成分進行二次處理，提出ICSA-PCA算法，在一定程度上解決了以上存在的問題。本文從空間統計學的角度，利用高光譜遙感影像的空間相關性，提出了一種定量的主成分選擇方法。

1 研究方法

空間統計學是建立在相鄰地理單元存在某種聯系的基本假設之上的統計學，將統計學和現代圖形計算技術結合起來，用直觀的方法展現空間數據中所隱含的空間分布、空間模式以及空間相互作用等特征[26]。地統計學是空間統計學的重要組成部分，而半變異函數理論是地統計學處理空間數據的方法，是探索展布于空間并呈現出一定的隨機性和結構性的自然現象的重要技術和方法，被用于描述空間相關性[27]。數據空間分布的相關性越大，即空間上聚集分布的現象越明顯。若所測值不表現出任何空間依賴關系，那么，這一變量表現出空間不相關性或空間隨機性。

變程、基臺值、塊金值是常用的表達空間自相關性的半變異函數基本參數。變程的大小反映了區域化變量影響范圍的大小，或者說反映該變量的自相關尺度。在變程距離之內，空間上越靠近在一起的點之間的相關性越大，相隔距離大于變程的點之間沒有自相關性。基臺值與塊金值之差表示由于采樣數據中存在空間自相關性引起的方差變化范圍，反映了數據的隨機性，稱為拱高。拱高和基臺值的比值可以反映數據空間相關性的強弱[28-29]。利用空間數據的以上特性可以很好地研究主成分選擇的問題。

本文利用變程、拱高/基臺值兩個半變異函數參數進行PCA后主成分的選擇，結合以上原理，基于空間統計學的高光譜遙感影像主成分選擇流程見圖1，主要過程如下：

1)對高光譜遙感影像數據集進行PCA變換，獲取主成分。在去除各主成分的趨勢項影響后對各主成分數據進行正態化轉換。

2)選定理論半變異函數模型對計算出的各主成分的實驗半變異函數進行擬合，從而獲取各半變異函數參數變程、基臺值、塊金值，由此計算出用于主成分選擇的變程、拱高/基臺值兩個參數。

3)聯合變程、拱高/基臺值進行主成分選擇，結合高光譜遙感影像實驗結果，從主觀和客觀兩個方面來綜合確定主成分選擇的經驗閾值。

4)利用支持向量機(support vector machine，SVM)算法進行分類，通過分類的Kappa系數、制圖精度、選擇的波段數等評價主成分選擇方法的效果[30-31]。與3種傳統選擇主成分的方法進行比較，評價本文提出的方法的好壞。

圖1 算法流程圖

2 實驗及其結果分析

2.1 仿真實驗

為便于研究，仿真數據大小設計為72×72。從美國地質調查局網站提供的地物波譜庫里隨機選擇4類地物波譜，按照規則格網設計仿真圖像(格網大小W=1,3,…,23)。為更加接近真實數據情況，添加信噪比分別為10，20，30，45的零均值高斯噪聲(圖2)。

圖2 仿真圖像(不同柵格大小W=1，3，…，23； 信噪比SNR=10，20，30，45)

如圖3、圖4所示，當W=1時，即相鄰的像元均為不同地物，圖像以隨機性為主，實驗半變異函數主要表現為塊金值。當W=3，5，…，23時，圖像的變程計算結果與設計的網格大小基本一致，反映了圖像空間相關性的范圍大小。拱高/基臺值計算結果受噪聲影響比較明顯，可以作為利用變程選擇主成分的輔助參數。仿真實驗從數據的空間相關性和隨機性兩個方面驗證了利用半變異函數參數拱高/基臺值、變程進行PCA后主成分選擇的有效性。

圖3 不同柵格大小仿真圖像計算的變程結果

圖4 不同柵格大小仿真圖像計算的拱高/基臺值結果

2.2 實際數據

2.2.1 實驗數據集

Indian Pines高光譜數據集是1992年由AVIRIS傳感器獲取的印第安納州西北部農業區的高光譜遙感影像數據的一部分，圖像大小為145×145像素，包含16類地物(圖5)。AVIRIS數據光譜范圍為0.4～2.45 μm，共224個波段，空間分辨率20 m。

(a) 影像(b) 真實標簽

ROSIS傳感器于2003年在意大利的北部Pavia大學獲取了2幅高光譜影像，University高光譜數據是該數據集的其中之一(圖6)。圖像大小為610×340 像素，空間分辨率1.3 m，包含9類地物。ROSIS傳感器共103個波段，光譜范圍為0.43～0.86 μm。

(a) 影像(b) 真實標簽

Salinas高光譜數據集由AVIRIS傳感器獲取的美國加利福尼亞州薩利納斯山谷區域。圖像大小為512×217像素，空間分辨率3.7 m，包含224個波段，該數據集地物類別包含16類(圖7)。

(a) 影像 (b) 真實標簽

研究中使用的Indian Pines，Pavia University(簡寫為Pavia U)和Salinas3種高光譜數據集獲取網站網址如下： http: //www.ehu.eus/ccwintco/index.php?title=Hyperspectral_Remote_Sensing_Scenes。

2.2.2 數據處理及結果

數據處理主要包含無信息波段的剔除、趨勢項去除、數據正態化、實驗半變異函數計算、理論半變異函數擬合等過程。

由于受水汽影響，Indian Pines數據集部分波段的成像效果比較差，本數據集中剔除的波段為104～108,150～163,220。同樣，Salina數據集剔除108～112,154～167,224波段。通過二階數據漂移估計，消除趨勢項影響。空間統計學中一般都假設數據是服從正態分布，本文在正態檢驗的基礎上，采用常態得分變換(normal score transform，NST)方法進行數據正態化的轉換，該方法相比傳統方法不受數據負值的影響[32]。

對于高光譜遙感影像的實驗半變異函數的計算，一般是分別計算圖像的0°,45°,90°和135°方向的實驗半變異函數曲線，研究其曲線變化特征，通過套和獲取最終的實驗半變異函數。選擇5種常用的模型——球狀模型、指數模型、高斯模型、線性有基臺值模型和普通線性模型模型，采用線性規劃法進行擬合，最終，選擇交叉驗證結果最佳的一個作為最終的理論模型[33-34]。一般采用交叉驗證后判定系數R2較大，或者殘差標準差較小的理論模型作為最終結果。實驗中采用判定系數R2和殘差標準差的比值作為理論模型的選擇標準。半變異函數參數計算結果如圖8所示，為便于圖面表達，圖示中舍棄了拱高/基臺值小于0.05的無意義主成分。

2.2.3 主成分選擇方法

從變程、拱高/基臺值的計算結果來看(圖8)，高光譜遙感影像PCA變換后編號較大的主成分主要表現為隨機噪聲，結果主要體現為塊金值，該特征與仿真實驗比較一致，主成分選擇中舍棄該類主成分。單獨利用變程或者拱高/基臺值也可以進行主成分的篩選，但是對于編號較大的主成分計算出的無意義結果不能很好的判斷。同時，變程和拱高/基臺值的結果具有明顯的互補性，對于一些無意義結果，同時通過2組參數可以有效的進行剔除。因此，本文提出綜合利用以上2組參數進行主成分選擇的思路。

2.2.4 閾值確定

利用變程、拱高/基臺值選擇主成分關鍵的問題就是閾值的確定。通過仿真實驗可以知道，圖像最小的空間相關性范圍大小為2，即相鄰像元是相關的，也就是變程為2。為便于研究，將變程增加到2.5作為對比實驗。拱高/基臺值體現了圖像的隨機性，一般認為當該值小于0.2～0.25時，數據表現為強的隨機性。為了便于研究，分別采用拱高/基臺值為0.2和0.25進行對比研究。在此基礎上，分別測試了變程=2、拱高/基臺值=0.2(表示為PC(2～0.2))； 變程=2、拱高/基臺值=0.25(表示為PC(2～0.25))； 變程=2.5、拱高/基臺值=0.2(表示為PC(2.5～0.2))； 變程=2.5、拱高/基臺值=0.25(表示為PC(2.5～0.25))的主成分選擇效果。

(a) Indian Pines數據集參數 (b) Pavia U數據集參數 (c) Salinas數據集參數

為了說明利用幾種閾值進行主成分選擇的效果，從主觀和客觀2個方面進行評價。主觀評價方法即觀察利用幾種閾值選擇出的主成分的圖像質量。以Indian Pines數據集為例，圖9為該數據集PCA變換后各主成分的縮略圖，表1為不同閾值篩選的主成分。圖9中排列順序為從左至右，從上至下，主成分編號依次增大，表示為PC1，PC2，…，PC200。實驗中發現，無論哪種閾值組合都可以剔除圖面質量較差的PC9。當拱高/基臺值固定時，增大變程會剔除更多的主成分。從表1結果來看，當拱高/基臺值=0.25時，PC103被剔除，同時，當變程由2增加到2.5會導致PC108，PC109被剔除。從圖9來看，PC103，PC108，PC109主成分的圖像細節仍然比較清楚。當變程固定時，減小拱高/基臺值會增加更多主要表現為隨機性的主成分。從表1結果來看，當變程=2時，拱高/基臺值由0.25減小到0.2會將PC22，PC29，PC55，PC59，PC103篩選進來。從圖9來看，PC22，PC29，PC59，PC103圖像細節仍然比較清楚，但是PC55圖像質量較差。因此還不能完全說明拱高/基臺值的閾值確定為哪個比較合適。

圖9 Indian Pines數據集PCA后各主成分縮略圖

表1 Indian Pines數據集不同閾值篩選的主成分

客觀評價方法是計算所選擇的主成分分類的Kappa系數進行對比評價。本文利用常用的SVM方法對幾種閾值的篩選結果進行分類，分類過程通過ENVI軟件實現，利用RBF核，設置gamma=0.1，penalty=100。利用幾組閾值選擇的主成分進行分類，從表2的分類精度結果來看均能得到較高的分類精度，從圖10的分類效果來看結果差別不明顯。綜合考慮主觀評價結果，利用PC(2～0.2)篩選出的主成分較多，利用PC(2.5～0.25)篩選出的主成分較少，PC(2.5～0.2)總體分類精度較高，因此，最終選擇變程=2.5，拱高/基臺值=0.2作為主成分選擇的閾值。

表2 不同閾值的Kappa系數

(a) PC(2～0.25)分類結果(b) PC(2～0.2)分類結果(c) PC(2.5～0.25)分類結果(d) PC(2.5～0.2)分類結果(e) 真實標簽

2.2.5 效果評價

為了進一步驗證提出的方法的有效性，采用了2種使用較為廣泛的傳統方法和一種創新方法進行對比研究。第一種是利用累積方差貢獻率進行主成分篩選，實驗中閾值定為0.99(表示為PC(0.99))[11]； 第二種是Rodarmel等[21]提出的可以用5%～10%的主成分個數進行分類研究，實驗中采用分類效果較好的10%(表示為PC(10%))； 第三種是Ibarrola-Ulzurrun等[22]創新提出的反映紋理特征指標即信息熵，選擇大于信息熵標準差的主成分進行分類(表示為PC(Entropy))。

各主成分選擇方法分類精度統計結果見表3，分類結果圖見圖11。對于Indian Pines和Salinas數據集，利用PC(2.5～0.2)分類的結果優于其他方法。對于Pavia U數據集，利用PC(2.5～0.2)分類的結果明顯優于PC(0.99)，且與PC(10%)、PC(Entropy)2種方法的分類精度差別不大。

表3 不同方法Kappa系數

(a) PC(0.99)分類結果(b) PC(10%)分類結果(c) PC(Entropy)分類結果(d) PC(2.5～0.2)分類結果(e) 真實標簽

表4為Indian Pines數據集各地類的制圖精度統計結果，數據為百分比。從表4可知，本文方法對于林地、大豆略耕地、燕麥地、牧草已割地、牧草地、玉米地、玉米未耕地和玉米略耕地的分類精度優于其他方法。同時，本文方法篩選出的主成分對于數量較少的地物較為敏感，塔樓、小麥地、燕麥地、牧草已割地、玉米地、苜蓿地、玉米地等分類效果明顯優于其他方法。表5為Pavia U數據集各地類的制圖精度統計結果，數據為百分比。從表5可知，本文方法對于裸地、柏油房頂、樹的分類精度優于其他方法。同時，本文方法篩選出的主成分對于數量較少的樹較為敏感，分類效果優于其他方法。表6為Salinas數據集各地類的制圖精度統計結果，數據為百分比。從表6可知，該數據集各地類總體分類精度比較高，本文方法對于未培育的葡萄園、長葉萵苣的分類精度優于其他方法。同時，本文方法篩選出的主成分對于類別較少的長葉萵苣_6wk地物較為敏感，分類效果優于其他方法。表7為不同方法所選擇的主成分個數。從表7可知，利用累計方差貢獻率大于0.99選擇的主成分個數因數據不同差別比較大，直接影響分類效果，不能作為一種適用于所有遙感影像數據的方法。利用10%的主成分個數的主成分開展分類效果可以，但無法解釋其物理意義，且結果受波段總數影響較大，結果具有隨機性。利用信息熵選擇的主成分個數因數據不同差別比較大，當地物易分類時，所選擇的主成分個數過多。利用信息熵進行主成分選擇會選擇出無信息的個別主成分，例如Indian Pines數據集的第49主成分，從圖9來看，該主成分無明顯的圖像信息。本文方法受數據影響較小，均能篩選出數量適中的主成分。本文方法不僅可以剔除一些編號雖然較小，但是圖像質量比較差的主成分，而且可以將編號較大，但是圖像質量較好的圖像參與運算。

表4 Indian Pines數據集制圖精度

表5 Pavia U數據集制圖精度

表6 Salinas數據集制圖精度

表7 不同數據集篩選出的主成分信息

3 結論

本文提出一種基于空間統計學的PCA變換后主成分選擇的新方法，利用半變異函數參數變程、拱高/基臺值的特性進行PCA變換后主成分的選擇，取得了理想的效果，得出如下結論：

1)仿真實驗證明了變程、拱高/基臺值可以有效表達高光譜遙感影像空間相關性的范圍和強弱。

2)變程、拱高/基臺值的結果具有明顯的互補性，聯合兩組參數可以有效剔除無意義主成分。

3)對比變程2～2.5和拱高/基臺值0.2～0.25的不同參數組合的分類結果，變程2.5、拱高/基臺值0.2的參數組合可以更加有效地篩選主成分。

4)和傳統方法相比，本文提出的方法可以剔除主成分編號較小，但是圖像質量較差的主成分，同時，篩選出主成分編號較大，但是圖像質量較好的主成分。

5)在基于分類的研究中，利用變程2.5、拱高/基臺值0.2進行PCA變換后主成分選擇，不僅能夠達到降維的目的，同時能夠保證足夠高的分類精度。和傳統方法相比，本文方法可以更好地識別數量比較少的地類。

由于半變異函數參數的計算也是一個研究比較多的問題，確定更加準確的半變異函數參數會對主成分的選擇產生一定的影響。除此之外，本文方法僅對PCA變換后的主成分選擇進行了探討，同時可以推廣到最大噪聲分數變換、獨立成分分析等降維方法中去。