高光譜圖像分類

1、機器學習課程項目報告高光譜圖像分類 基于CNN和ELM 學 院 信息工程學院 專 業(yè) 電子與通信工程 學 號 2111603035 學生姓名 曹發(fā)賢 同組學生 陳惠明、陳濤 碩士導師 楊志景 2016 年 11 月一、 項目意義與價值高光譜遙感技術(shù)起源于 20 世紀 80年代初，是在多光譜遙感技術(shù)基礎(chǔ)之上發(fā)展起來的1。高光譜遙感能夠通過成像光譜儀在可見光、近紅外、短波紅外、中紅外等電磁波譜范圍獲取近似連續(xù)的光譜曲線，將表征地物幾何位置關(guān)系的空間信息與表征地物屬性特征的光譜信息有機地融合在了一起，使得提取地物的細節(jié)信息成為可能。隨著新型成像光譜儀的光譜分辨率的提高，人們對相關(guān)地物的光譜屬性特征的

2、了解也不斷深入，許多隱藏在狹窄光譜范圍內(nèi)的地物特性逐漸被人們所發(fā)現(xiàn)，這些因素大大加速了遙感技術(shù)的發(fā)展，使高光譜遙感成為 21 世紀遙感技術(shù)領(lǐng)域重要的研究方向之一。在將高光譜數(shù)據(jù)應用于各領(lǐng)域之前，必須進行必要的數(shù)據(jù)處理。常用的數(shù)據(jù)處理技術(shù)方法包括：數(shù)據(jù)降維、目標檢測、變化檢測等。其中，分類是遙感數(shù)據(jù)處理中比較重要的環(huán)節(jié)，分類結(jié)果不但直接提取了影像數(shù)據(jù)有效信息，可以直接運用于實際需求中，同時也是實現(xiàn)各種應用的前提，為后續(xù)應用提供有用的數(shù)據(jù)信息和技術(shù)支持，如為目標檢測提供先驗信息、為解混合提供端元信息等。相對于多光譜遙感而言，由于高光譜遙感的波譜覆蓋范圍較寬，因此我們可以根據(jù)需要選擇特定的波段來突顯

3、地物特征，從而能夠精確地處理地物的光譜信2。目前，許多國家開展大量的科研項目對高光譜遙感進行研究，研制出許多不同類型的成像光譜儀。高光譜遙感正逐步從地面遙感發(fā)展到航空遙感和航天遙感，并在地圖繪制、資源勘探、農(nóng)作物監(jiān)測、精細農(nóng)業(yè)、海洋環(huán)境監(jiān)測等領(lǐng)域發(fā)揮重要的作用。高光譜遙感技術(shù)雖然是遙感領(lǐng)域的新技術(shù)，但是高光譜圖像的分類一直制約著高光譜遙感的應用3,4，因此對其進行研究顯得尤為重要。高光譜遙感圖像較高的光譜分辨率給傳統(tǒng)的圖像分類識別算法提出嚴峻的挑戰(zhàn)。波段維數(shù)的增加不僅加重了數(shù)據(jù)的存儲與傳輸?shù)呢摀?，同時也加劇了數(shù)據(jù)處理過程的復雜性，并且由于波段與波段間存在著大量的冗余信息，從而使得傳統(tǒng)圖像分類算

4、法并不適用于高光譜遙感圖像的分類。傳統(tǒng)的分類方法往往需要很多數(shù)目的已知類別的訓練樣本，從而導致計算量大，時間復雜度高。另外，如果訓練樣本的數(shù)目較少，采用傳統(tǒng)分類算法進行分類時分類精度往往是很低的，因此提高分類精度并減少運算量成為高光譜領(lǐng)域的熱點問題。高光譜遙感圖像的波段數(shù)目多，并且波段與波段間存在著很大的相關(guān)性，因此在進行遙感圖像的處理（例如分類）之前通常需要進行降維預處理，這樣做不僅可以降低數(shù)據(jù)空間的維數(shù)，減少冗余信息，而且還有利于人工圖像解譯和后續(xù)分類處理和地物識別，從而為解決高光譜遙感分類的難點提供了方便5。二、 高光譜圖像分類的發(fā)展與現(xiàn)狀高光譜圖像分類作為高光譜圖像的基礎(chǔ)研究，一直是高

5、光譜圖像重要的信息獲取手段，它的主要目標是根據(jù)待測地物的空間幾何信息與光譜信息將圖像中的每個像素劃分為不同的類別。高光譜圖像分類按照是否有已知類別的訓練樣本的參與，高光譜圖像的分類方式分為監(jiān)督分類與非監(jiān)督分類6。在遙感圖像自動分類中，傳統(tǒng)的基于數(shù)理統(tǒng)計的分類方法，主要包括最小距離分類、最大似然分類、波譜角分類、混合距離法分類等，主要依賴地物的光譜屬性，基于單個像元進行分類。統(tǒng)計模式識別方法本身的不足：1、最大似然法計算強度大，且要求數(shù)據(jù)服從正態(tài)分布2、K-means聚類分類精度低，分類精度依賴于初始聚類中心3、最小距離法沒有考慮各類別的協(xié)方差矩陣，對訓練樣本數(shù)目要求低近年來對于神經(jīng)網(wǎng)絡(luò)分類方法

6、的研究相當活躍。它區(qū)別于傳統(tǒng)的分類方法：在處理模式分類問題時，并不基于某個假定的概率分布，在無監(jiān)督分類中，從特征空間到模式空間的映射是通過網(wǎng)絡(luò)自組織完成的。在監(jiān)督分類中，網(wǎng)絡(luò)通過對訓練樣本的學習，獲得權(quán)值，形成分類器，且具備容錯性。人工神經(jīng)網(wǎng)絡(luò)（ANN）分類方法一般可以獲得更高精度的分類結(jié)果，因此ANN方法在遙感分類中被廣泛應用，特別是對于復雜類型的地物類型分類，ANN方法顯示了其優(yōu)越性。專家系統(tǒng)分類法也在遙感分類取得了一定的應用。專家系統(tǒng)是模擬人類邏輯思維的智能系統(tǒng)，將其應用于遙感分類最大的優(yōu)點就是可以充分利用更多的輔助分類數(shù)據(jù)。不過由于專家系統(tǒng)知識庫難以建立，影響了它的進一步發(fā)展。支持向量

7、機（SVM）具有嚴格的理論基礎(chǔ)，能較好地解決小樣本、非線性、高維數(shù)等問題，被成功地應用到多光譜、高光譜遙感圖像分類領(lǐng)域。對于高光譜數(shù)據(jù)而言，由于波段多、數(shù)據(jù)量大、數(shù)據(jù)不確定性等，易受Hughes現(xiàn)象（即訓練樣本固定時，分類精度隨特征維數(shù)的增加而下降）影響。而樣本的獲取在高光譜分類中往往是一項比較困難的工作，特別是采用高維特征向量時要求每類的樣本數(shù)都要比特征維數(shù)高，因此在高維信息處理中的精度與效率和高光譜遙感信息精細光譜與大數(shù)據(jù)量之間仍然存在著極大的矛盾。三、 卷積神經(jīng)網(wǎng)絡(luò)理論基礎(chǔ)卷積神經(jīng)網(wǎng)絡(luò)是人工神經(jīng)網(wǎng)絡(luò)的一種，它的權(quán)值共享網(wǎng)絡(luò)結(jié)構(gòu)使之更類似于生物神經(jīng)網(wǎng)絡(luò)，降低了網(wǎng)絡(luò)模型的復雜度，減少的權(quán)值的

8、數(shù)量以節(jié)約訓練和測試的計算時間。該優(yōu)點在網(wǎng)絡(luò)的輸入是多維圖像時表現(xiàn)得更加明顯，使圖像可以直接作為網(wǎng)絡(luò)的輸入，避免了傳統(tǒng)識別算法中復雜的特征提取和數(shù)據(jù)重建過程。卷積神經(jīng)網(wǎng)絡(luò)是為識別二維數(shù)據(jù)而專門設(shè)計的一個多層感知機，其網(wǎng)絡(luò)對平移、比例變化和傾斜等具有高度不變性7。在CNN中，圖像的一小部分（局部感受區(qū)域）作為層級結(jié)構(gòu)的最低層的輸入，信息再依次傳輸?shù)讲煌膶?，每層通過一個數(shù)字濾波器去獲得觀測數(shù)據(jù)的最顯著的特征。這個方法能夠獲取對平移、縮放和旋轉(zhuǎn)不變的觀測數(shù)據(jù)的顯著特征，因為圖像的局部感受區(qū)域允許神經(jīng)元或者處理單元可以訪問到最基礎(chǔ)的特征，例如定向邊緣或者角點。卷積神經(jīng)網(wǎng)絡(luò)是一個多層的神經(jīng)網(wǎng)絡(luò)，每層由

9、多個二維平面組成，每個平面又由多個獨立的神經(jīng)元組成。圖2為卷積神經(jīng)網(wǎng)絡(luò)的整體結(jié)構(gòu)圖。一般地，C層（卷積層）為特征提取層，每個神經(jīng)元的輸入與前一層的局部感受野相連，并提取該局部的特征，一旦該局部特征被提取后，它與其他特征間的位置關(guān)系也隨之確定下來；S層（下采樣層）是特征映射層，網(wǎng)絡(luò)的每個計算層由多個特征映射組成，每個特征映射為一個平面，平面上所有神經(jīng)元的權(quán)值相等。特征映射結(jié)構(gòu)采用sigmoid函數(shù)等映射函數(shù)作為卷積網(wǎng)絡(luò)的激活函數(shù)，使得特征映射具有位移不變性。此外，由于一個映射面上的神經(jīng)元共享權(quán)值，因而減少了網(wǎng)絡(luò)自由參數(shù)的個數(shù)，降低了網(wǎng)絡(luò)參數(shù)選擇的復雜度。卷積神經(jīng)網(wǎng)絡(luò)中的每一個特征提取層（C-層）

10、都緊跟著一個用來求局部平均與二次提取的計算層（S-層），這種特有的兩次特征提取結(jié)構(gòu)使網(wǎng)絡(luò)在識別時對輸入樣本有較高的畸變?nèi)萑棠芰?。卷積神經(jīng)網(wǎng)絡(luò)采用有監(jiān)督學習的方式進行訓練，即任何一個訓練樣本的類別是已知的，訓練樣本在空間中的分布不再是依據(jù)其自然分布傾向來劃分，而是根據(jù)同類樣本和不同類樣本中的空間關(guān)系進行劃分，這需要不斷調(diào)整網(wǎng)絡(luò)模型的參數(shù)用以劃分樣本空間的分類邊界的位置，是一個耗時且復雜的學習訓練過程8。神經(jīng)網(wǎng)絡(luò)在進行訓練時，所有的網(wǎng)絡(luò)權(quán)值都用一些不同的小隨機數(shù)進行初始化，這些小的隨機數(shù)能偶保證網(wǎng)絡(luò)不會因為權(quán)值過大而進入飽和狀態(tài)，導致訓練失敗。神經(jīng)網(wǎng)絡(luò)訓練算法包括4個主要部分：（1）樣本集中取出樣

11、本 (X,yp) 并將其輸入網(wǎng)絡(luò)，X代表圖像數(shù)組，yp代表其類別；（2）計算此次輸入相應的實際輸出 Op ,這是一個前向傳播的過程；（3）用一個指定的損失函數(shù)計算出實際輸出 Op 與理想輸出 Yp 的誤差；（4）按極小化誤差的方法反向傳播調(diào)整網(wǎng)絡(luò)權(quán)值。四、 極限學習機極限學習機（extreme learning machine）ELM是一種簡單易用、有效的單隱層前饋神經(jīng)網(wǎng)絡(luò)SLFNs學習算法。2004年由南洋理工大學黃廣斌副教授提出。傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)學習算法（如BP算法）需要人為設(shè)置大量的網(wǎng)絡(luò)訓練參數(shù)，并且很容易產(chǎn)生局部最優(yōu)解。極限學習機只需要設(shè)置網(wǎng)絡(luò)的隱層節(jié)點個數(shù)，在算法執(zhí)行過程中不需要調(diào)整網(wǎng)

12、絡(luò)的輸入權(quán)值以及隱元的偏置，并且產(chǎn)生唯一的最優(yōu)解，因此具有學習速度快且泛化性能好的優(yōu)點。極限學習機的網(wǎng)絡(luò)訓練模型采用前向單隱層結(jié)構(gòu)。設(shè) 分別為網(wǎng)絡(luò)輸入層、隱含層和輸出層的節(jié)點數(shù)， 是隱層神經(jīng)元的激活函數(shù)， 為閾值。設(shè)有 個不同樣本 ， ，其中 ，則極限學習機的網(wǎng)絡(luò)訓練模型如圖1所示。圖1 極限學習機的網(wǎng)絡(luò)訓練模型極限學習機的網(wǎng)絡(luò)模型可用數(shù)學表達式表示如下：式中， 表示連接網(wǎng)絡(luò)輸入層節(jié)點與第i個隱層節(jié)點的輸入權(quán)值向量； 表示連接第i個隱層節(jié)點與網(wǎng)絡(luò)輸出層節(jié)點的輸出權(quán)值向量； 表示網(wǎng)絡(luò)輸出值。極限學習機的代價函數(shù)E可表示為式中， ，包含了網(wǎng)絡(luò)輸入權(quán)值及隱層節(jié)點閾值。Huang等指出極限學習機的懸鏈

13、目標就是尋求最優(yōu)的S，使得網(wǎng)絡(luò)輸出值與對應實際值誤差最小，即 ?？蛇M一步寫為式中，H表示網(wǎng)絡(luò)關(guān)于樣本的隱層輸出矩陣，表示輸出權(quán)值矩陣，T表示樣本集的目標值矩陣，H，T分別定義如下：極限學習機的網(wǎng)絡(luò)訓練過程可歸結(jié)為一個非線性優(yōu)化問題。當網(wǎng)絡(luò)隱層節(jié)點的激活函數(shù)無限可微時，網(wǎng)絡(luò)的輸入權(quán)值和隱層節(jié)點閾值可隨機賦值，此時矩陣H為一常數(shù)矩陣，極限學習機的學習過程可等價為求取線性系統(tǒng) 最小范數(shù)的最小二乘解 ，其計算式為式中時矩陣的MP廣義逆。五、 具體實現(xiàn)及主要代碼1、 訓練的樣本及其樣本圖 2、原圖及進行分類后的圖和各個高光譜數(shù)據(jù)集每一類的分類精度The result of Pavia Universit

14、y The result of Pavia Center The result of Salinas 3、主要代碼The main code of Pavia University of Matlab% Convolution Neural Network And Extrem Learning Machine%Test 1 for PaviaUniversity% STEP 1: InitializecnnConfig=config();theta, meta = cnnInitParams(cnnConfig); n_class = cnnConfig.layermeta.numLayer

15、s.dimension;load PaviaU;load PaviaU_gt;images, labels, indexs,samimage = loadtrain(paviaU, paviaU_gt);d = cnnConfig.layer1.dimension;images = reshape(images,d(1),d(2),d(3),);% STEP 2: Train (The Cnn And Elm) options.epochs = 1;options.minibatch = 30;options.alpha = 0.05;newtheta = minFuncSGD(x,y,z)

16、cnnCost(x,y,z,cnnConfig,meta),theta,images,labels,options);K = cnnExtract(newtheta,images,cnnConfig,meta);TrainingTime,TrainingAccuracy,InputWeight,BiasofHiddenNeurons,OutputWeight,NumberofOutputNeurons =elmtrain(K,labels ,1, 900);% STEP 3: TesttestImages, testLabels, testIndexs = loadtest(paviaU, p

17、aviaU_gt);testImages = reshape(testImages,d(1),d(2),d(3),);row, col = size(paviaU_gt);testK = cnnExtract(newtheta,testImages,cnnConfig,meta);TestingTime, TestingAccuracy,testoutputlabel,actualoutputs = elmpredict(testK,testLabels,testIndexs,1,InputWeight,BiasofHiddenNeurons,OutputWeight,NumberofOutp

18、utNeurons,row,col);predimage1=zeros(row,col);predimage1(testIndexs)=testoutputlabel;figure,imagesc(predimage1);axis off;OA,kappa,AA,CA= calcError(testLabels, predimage1(testIndexs), 1:n_class);fprintf(cnn+elm Overall Accuracy is %fn,OA);fprintf(cnn+elm Average Accuracy is %fn,AA);fprintf(cnn+elm CA

19、Acuuyracy is %fn,CA);fprintf(cnn+elm kappa is %fn,kappa);figure,imagesc(paviaU_gt);axis off;The main code of Pavia Center of Matlab% Convolution Neural Network And Extrem Learning Machine%Test 2 for Pavia centercnnConfig = config3();theta, meta = cnnInitParams(cnnConfig);n_class = cnnConfig.layermet

20、a.numLayers.dimension;load Pavia;load Pavia_gt;images, labels, indexs = loadtrain3(pavia,pavia_gt);d = cnnConfig.layer1.dimension;images = reshape(images,d(1),d(2),d(3),);% STEP 2: Train (The Cnn And Elm) options.epochs = 1;options.minibatch = 25;options.alpha = 0.05;newtheta = minFuncSGD(x,y,z) cnn

21、Cost(x,y,z,cnnConfig,meta),theta,images,labels,options);K = cnnExtract(newtheta,images,cnnConfig,meta);TrainingTime,TrainingAccuracy,InputWeight,BiasofHiddenNeurons,OutputWeight,NumberofOutputNeurons =elmtrain(K,labels ,1, 850);% STEP 3: TesttestImages, testLabels, testIndexs = loadtest3(pavia, pavi

22、a_gt);testImages = reshape(testImages,d(1),d(2),d(3),);row, col = size(pavia_gt);testK = cnnExtract(newtheta,testImages,cnnConfig,meta);TestingTime, TestingAccuracy,testoutputlabel = elmpredict(testK,testLabels,1,InputWeight,BiasofHiddenNeurons,OutputWeight,NumberofOutputNeurons,row,col);predimage1=

23、zeros(row,col);predimage1(testIndexs)=testoutputlabel;figure,imagesc(predimage1);axis offOA,kappa,AA,CA= calcError(testLabels, predimage1(testIndexs), 1:n_class);fprintf(cnn+elm Overall Accuracy is %fn,OA);fprintf(cnn+elm Average Accuracy is %fn,AA);fprintf(cnn+elm CA Acuuyracy is %fn,CA);fprintf(cn

24、n+elm kappa is %fn,kappa);figure,imagesc(pavia_gt);axis offThe main code of Salinas of Matlab% Convolution Neural Network And Extrem Learning Machine%Test 1 for Salinas% STEP 1: InitializecnnConfig=config();theta, meta = cnnInitParams(cnnConfig);n_class = cnnConfig.layermeta.numLayers.dimension;load

25、 Salinas_corrected;load Salinas_gt;images, labels, indexs,samimage = loadtrain(salinas_corrected, salinas_gt);d = cnnConfig.layer1.dimension;images = reshape(images,d(1),d(2),d(3),);% STEP 2: Train (The Cnn And Elm) options.epochs = 1;options.minibatch = 30;options.alpha = 0.05;newtheta = minFuncSGD

26、(x,y,z) cnnCost(x,y,z,cnnConfig,meta),theta,images,labels,options);K = cnnExtract(newtheta,images,cnnConfig,meta);TrainingTime,TrainingAccuracy,InputWeight,BiasofHiddenNeurons,OutputWeight,NumberofOutputNeurons =elmtrain(K,labels ,1, 2300);% STEP 3: TesttestImages, testLabels, testIndexs = loadtest(

27、salinas_corrected,salinas_gt);testImages = reshape(testImages,d(1),d(2),d(3),);row, col = size(salinas_gt);testK = cnnExtract(newtheta,testImages,cnnConfig,meta);TestingTime, TestingAccuracy,testoutputlabel,actualoutputs = elmpredict(testK,testLabels,testIndexs,1,InputWeight,BiasofHiddenNeurons,Outp

28、utWeight,NumberofOutputNeurons,row,col);predimage1=zeros(row,col);predimage1(testIndexs)=testoutputlabel;figure,imagesc(predimage1);axis off;OA,kappa,AA,CA= calcError(testLabels, predimage1(testIndexs), 1:n_class);fprintf(cnn+elm Overall Accuracy is %fn,OA);fprintf(cnn+elm Average Accuracy is %fn,AA

29、);fprintf(cnn+elm CA Acuuyracy is %fn,CA);fprintf(cnn+elm kappa is %fn,kappa);figure,imagesc(salinas_gt);axis off;六、 參考文獻1 童慶禧，張兵，鄭蘭芬.高光譜遙感原理.技術(shù)與應用M.北京：高等教育出版社, 2006 2 S.Wang,C. I. Chang.Variable-Number Variable-Band Selection for Feature Characterization in Hyperspectral Signatures J.IEEE Transanctions on Geoscience and Remote Sensing, 2007, 45(9): 2979-2992 3 Q. Du, C. I. Chang.Linear mixture analysis-based compression for hyperspectral