高光谱图像分类

《机器学习》课程项目报告
高光谱图像分类
——基于CNN和ELM
学院信息工程学院
专业电子与通信工程
学号35
学生姓名曹发贤
同组学生陈惠明、陈涛
硕士导师 _______ 杨志景_______
2016年11月
一、项目意义与价值
高光谱遥感技术起源于20世纪80年代初，是在多光谱遥感技术基础之上发展起来的［11 o高光谱遥感能够通过成像光谱仪在可见光、近红外、短波红外、中红外等电磁波谱范围获取近似连续的光谱曲线，将表征地物儿何位置关系的空间信息与表征地物属性特征的光谱信息有机地融合在了一起，使得提取地物的细节信息成为可能。

随着新型成像光谱仪的光谱分辨率的提高，人们对相关地物的光谱属性特征的了解也不断深入，许多隐藏在狭窄光谱范围内的地物特性逐渐被人们所发现，这些因素大大加速了遥感技术的发展，使高光谱遥感成为21 世纪遥感技术领域重要的研究方向之一。

在将高光谱数据应用于各领域之前，必须进行必要的数据处理。

常用的数据处理技术方法包括：数据降维、目标检测、变化检测等。

其中，分类是遥感数据处理中比较重要的环节，分类结果不但直接提取了影像数据有效信息，可以直接运用于实际需求中，同时也是实现各种应用的前提，为后续应用提供有用的数据信息和技术支持，如为目标检测提供先验信息、为解混合提供端元信息等。

相对于多光谱遥感而言，由于高光谱遥感的波谱覆盖范围较宽，因此我们可以根据需要选择特定的波段来突显地物特征，从而能够精确地处理地物的光谱信［2］。

目前，许多国家开展大量的科研项目对高光谱遥感进行研究，研制出许多不同类型的成像光谱仪。

高光谱遥感正逐步从地面遥感发展到航空遥感和航天遥感，并在地图绘制、资源勘
探、农作物监测、精细农业、海洋环境监测等领域发挥重要的作用。

高光谱遥感技术虽然是遥感领域的新技术，但是高光谱图像的分类一直制约着高光谱遥感的应用［3, 4］,因此对其进行研究显得尤为重要。

高光谱遥感图像较高的光谱分辨率给传统的图像分类识别算法提出严峻的挑战。

波段维数的增加不仅加重了数据的存储与传输的负担, 同时也加剧了数据处理过程的复杂性，并且由于波段与波段间存在着大量的兀余信息，从而使得传统图像分类算法并不适用于高光谱遥感图像的分类。

传统的分类方法往往需要很多数目的己知类别的训练样本，从而导致计算量大，时间复杂度高。

另外，如果训练样本的数目较少，采用传统分类算法进行分类时分类精度往往是很低的，因此提高分类精度并减少运算量成为高光谱领域的热点问题。

高光谱遥感图像的波段数目多，并且波段与波段间存在着很大的相关性，因此在进行遥感图像的处理（例如分类）之前通常需要进行降维预处理，这样做不仅可以降低数据空间的维数，减少冗余信息，而且还有利于人工图像解译和后续分类处理和地物识别，从而为解决高光谱遥感分类的难点提供了方便［5］。

二、高光谱图像分类的发展与现状
高光谱图像分类作为高光谱图像的基础研究，一直是高光谱图像重要的信息获取手段，它的主要目标是根据待测地物的空间儿何信息与光谱信息将图像中的每个像素划分为不同的类别。

高光谱图像分类按照是否有已
知类别的训练样本的参与，高光谱图像的分类方式分为监督分类与非监督分类［6］。

在遥感图像自动分类中，传统的基于数理统计的分类方法，主要包括最小距离分类、最大似然分类、波谱角分类、混合距离法分类等，主要依赖地物的光谱属性，基于单个像元进行分类。

统计模式识别方法本身的不足：
1、最大似然法计算强度大，且要求数据服从正态分布
2、K-means聚类分类精度低，分类精度依赖于初始聚类中心
3、最小距离法没有考虑各类别的协方差矩阵，对训练样本数目要求低
近年来对于神经网络分类方法的研究相当活跃。

它区别于传统的分类方法：在处理模式分类问题时，并不基于某个假定的概率分布，在无监督分类中，从特征空间到模式空间的映射是通过网络自组织完成的。

在监督分类中，网络通过对训练样本的学习，获得权值，形成分类器，且具备容错性。

人工神经网络（ANN）分类方法一般可以获得更高精度的分类结果，因此ANN方法在遥感分类中被广泛应用，特别是对于复杂类型的地物类型分类，ANN方法显示了其优越性。

专家系统分类法也在遥感分类取得了一定的应用。

专家系统是模拟人类逻辑思维的智能系统，将其应用于遥感分类最大的优点就是可以充分利用更多的辅助分类数据。

不过由于专家系统知识库难以建立，影响了它的进一步发展。

支持向量机（SVM）具有严格的理论基础，能较好地解决小样本、非对于
高光谱数据而言，由于波段多、数据量大、数据不确定性等，易受Hughes现象（即训练样本固定时，分类精度随特征维数的增加而下降）影响。

而样本的获取在高光谱分类中往往是一项比较困难的工作，特别是采用高维特征向量时要求每类的样本数都要比特征维数高，因此在高维信息处理中的精度与效率和高光谱遥感信息精细光谱与大数据量之间仍然存在着极大的矛盾。

线性、高维数等问题，被成功地应用到多光谱、高光谱遥感图像分类领域。

个平面又由多个独立的神经元组成。

图2为卷积神经网络的整体结构图。

三、卷积神经网络理论基础
卷积神经网络是人工神经网络的一种，它的权值共享网络结构使之更 类似
于生物神经网络，降低了网络模型的复杂度，减少的权值的数量以节 约训练和测试的计算时间。

该优点在网络的输入是多维图像时表现得更加 明显，使图像可以直接作为网络的输入，避免了传统识别算法中复杂的特 征提取和数据重建过程。

卷积神经网络是为识别二维数据而专门设计的一 个多层感知机，其网络对平移、比例变化和倾斜等具有高度不变性［7］。

在CNN 中，图像的一小部分（局部感受区域）作为层级结构的最低层 的输
入，信息再依次传输到不同的层，每层通过一个数字滤波器去获得观 测数据的最显着的特征。

这个方法能够获取对平移、缩放和旋转不变的观 测数据的显着特征，因为图像的局部感受区域允许神经元或者处理单元可 以访问到最基础的特征，例如定向边缘或者角点。

卷积神经网络是一个多层的神经网络，每层由多个二维平面组成，每
一般地，C 层（卷积层）为特征提取层，每个神经元的输入与前一层的 局
部感受野相连，并提取该局部的特征，一旦该局部特征被提取后，它与
其他特
征间的位置关系也随之确定下来；s层（下采样层）是特征映射层，网络的每个计算层由多个特征映射组成，每个特征映射为一个平面，平面上所有神经元的权值相等。

特征映射结构釆用sigmoid函数等映射函数作为卷积网络的激活函数，使得特征映射具有位移不变性。

此外，由于一个映射面上的神经元共享权值，因而减少了网络自由参数的个数，降低了网络参数选择的复杂度。

卷积神经网络中的每一个特征提取层（C-层）都紧跟着一个用來求局部平均与二次提取的计算层（S-层）, 这种特有的两次特征提取结构使网络在识別时对输入样本有较高的畸变容忍能力。

卷积神经网络采用有监督学习的方式进行训练，即任何一个训练样本的类别是已知的，训练样本在空间中的分布不再是依据其自然分布倾向來
划分，而是根据同类样本和不同类样本中的空间关系进行划分，这需要不
断调整网络模型的参数用以划分样本空间的分类边界的位置，是一个耗时且复杂的学习训练过程［8］。

神经网络在进行训练时，所有的网络权值都用一些不同的小随机数进行初始化，这些小的随机数能偶保证网络不会因为权值过大而进入饱和状态，导致训练失败。

神经网络训练算法包括4个主要部分：
（1）样本集中取出样本（X,*）并将其输入网络，X代表图像数组，yp 代表其类别；
（2）计算此次输入相应的实际输出Op，这是一个前向传播的过程；
（3）用一个指定的损失函数计算出实际输出Op与理想输出Yp的误差；
（4）按极小化误差的方法反向传播调整网络权值。

四、极限学习机
极限学习机（extreme learning machine） ELM是一种简单易用、有效的单隐层前馈神经网络SLFNs学习算法。

2004年由南洋理工大学黄广斌副教授提出。

传统的神经网络学习算法（如BP算法）需要人为设置大量的网络训练参数，并且很容易产生局部最优解。

极限学习机只需要设置网络的隐层节点个数，在算法执行过程中不需要调整网络的输入权值以及隐元的偏置，并且产生唯一的最优解，因此具有学习速度快且泛化性能好的优点。

极限学习机的网络训练模型采用前向单隐层结构。

设uM,”分别为网络输入层、隐含层和输出层的节点数，g（x）是隐层神经元的激活函数，
b,为阈值。

设有N 个不同样本(%,,/,) , \5iSN ,其中
兀=［无i,兀2,x in,丫 e R"\t, = ［t a,切,•••，订已R"，则极限学习机的网络训练模型如图1所示。

图1极限学习机的网络训练模型
极限学习机的网络模型可用数学表达式表示如下：
M
X 卩ig㈣吒+b) = Oj,j =1,2,…,N
式中，表示连接网络输入层节点与第八个隐层节点的输入权值向量；0尸\虫、际…朋表示连接第：个隐层节点与网络输出层节点的输出权值向量；=［。

门,o i2,『表示网络输出值。

极限学习机的代价函数疋可表示为
N
E(S，0) =丈加-引
式中，s = (q,Q,i = l,2,...,M)，包含了网络输入权值及隐层节点阈值o Huang 等指出极限学习机的悬链目标就是寻求最优的S, 0,使得网络输出值与对应实际值误差最小，即min(E(S,0))。

min(E(S,0))可进一步写为
式中，〃表示网络关于样本的隐层输出矩阵，0表示输出权值矩阵，T 表 示样本集的目标値矩阵，H. 0, 7分别定义如下:
g （QZ+勺）…g （®a+如）
g （©心+勺）…&（©血+如）
极限学习机的网络训练过程可归结为一个非线性优化问题。

当网络隐 层节点的激活函数无限可微时，网络的输入权值和隐层节点阈值可随机赋 值，此时矩阵H 为一常数矩阵，极限学习机的学习过程可等价为求取线性 系统Hp = T 最小范数的最小二乘解B ，其计算式为
B = H*T
式中H+时矩阵H 的MP 广义逆。

五、具体实现及主要代码
1、 训练的样本及其样本图
min （E （S,0）） =
min ||//（①，…®…厶f 心…旳）0-口
ROSIS Pavia University
No Name Train Test
1 J Asphalt—5486631
2 n Meadows64018649
3 J Gravel392 _J2099
4Trees3064
5 J Metal sheets2帘—1345
6 J Bare soil532 I5029
7 J Bitumen3751330
8 J Bricks5143682
9Shadows231 ___ 947
Total392142776
ROSIS Pavia Center
No Name Train Test
Water82465971 2Trees8207598
3Meadow8243090
4Bricks8082685
5Soil8206584
6Asphalt8169284
7Bitumen8087287
8Tile126042826
9Shadows4762863
Total7456148152
Salinas
No Name Train Test
1 n Brocoli_green_weed_l3012009
2 J Brocoli_green_weed_25593726 |
3Fallow2961976
4 J Fallow_rough_plow209
5Fallow_smooth4022678
6 J Stubble5943959
7Celery5373579
8 J Grapes^u ntrained1691辻27久
9 n Soil_vinyard_develop9306203
10Corn_sensced_green_weeds4923278
1 11Lettuce」。

maine_4wk1601068
12Lettuce」o m a i ne_5w k2891927
1 13Lettu ce_r o m a i ne_6w k137916
14Lettuce」o m a i ne_7w k16i1070
15Vinyard_untrai ned10907268
16Vinyard.verticaLtrellis271T807
Total811954129
羽dj b wma宴o目IQ J己心朋$ u也i0
50 100 150 2(0 ^0 3X1 100 20) 300 400 500 9)0
700
file Edit V»tw Insert Tools Cfeskicp Wrdow Help
□ J k・・m・龙D0 H0
2、原图及进行分类后的图和各个高光谱数据集每一类的分类精度
The result of Pavia University
fife Edit
\伽 矗ert Took Cesidoo Xtido*
Hdp
Odd ；
k 「祷旅存
Re Edit Viw xsert Took Desktop Wrriw Hdp
Iestinglme -
1.6536
hsti^Aozuraoy =
刪
cnn+eln Overall Accura^ is 0.932392 cnn+eln Averse Accura^ is 0. 939845 cnn+eh CA Aaxyracy is 0.S92626 cnn+eln CA Aoaiyracy is 0.942946 cnn+eh CA A'^nyracy is 0.S64697 cnn+eln CA Aajuyracy is 0.980418 cim+eln CA Aajuyracy is 0.987361 cnn+eh CA Acwqrrat/ is 0.927620 cnn+eln CA AaiTiacy is 0.954135
cnn+eln CA AajiTracy is 0.910918 cnn+eh CA Aaiqrracy is 0.&97S88 cnn+eln kappa is 0.911218
Training Time =
1. 9344
TrainingAccuracy =
O ・ 9707
第鹃b X ◎驛盒Q 0@ iQ
The result of Pavia Center
lest inglime =
5.5692
TestingAccuracy =
0.9892
IrainingTime =
3. 3384 TrainingAccuracy =
0. 9862cnn+ela Overall Accuracy is 0.989227 cnn+ela
cnn+ela
cnn+el>
cnn+elm
cnn+elB
cnn+ela
cnn+ela
cnn+ela
cnn+ela
cnn+ela
Average Accuracy is 0. 976868
CA
CA
CA
CA
CA
CA
CA
CA
CA
Acuuyracy
Acuuyracy
Acuuyracy
Acuuyracy
Acuuyracy
Acuuyracy
Acuuyracy
Acuuyracy
Acuuyracy
0.999803
0.958147
0.970550
0.978399
is 0.965067
0.982483
0.946068
0.992341
0.998952
cnn+ela kappa is 0.984758
The resuIt of Salinas
5.6316
Iesting Ac curacy =
0.9490
3、主要代码
The main code of Pavia University of Matlab
%% Convolution Neural Network And Extrem Learning Machine
Irauungline =
25.4282
Iranung Accuracy
=
0.9fiSi
cnn+eln.
CA AcoiiyrGcy
is 0. 999502 cnn+eln CA Acuuyracy as 0. 998121 cnn+eln CA Acuuyr2scy is 0. 998988 cnn+eln. CA Acuuyr2scy 2X 0. 997131 cnn+eln. CA Acuuyracy 13 0. 995892 cnn+eln. CA Acuuyracy 13 0. 962112 cnn+elD. CA Acuuyracy 13 0. 996927 crm+cln. CA Acuuyracy 13 0. 904001 crm+eln. CA Acuuyracy is 0. 996292 cnn+eln. CA Acuuyracy is 0. 966138 cnn+eln. CA Acuuyracy is 0.993446 cnn+eln. CA Acmiyrscy is I. 000000 cnn+ela CA Acuuyrscy is 0. 993450 cnn+eln CA AcoiiyrGcy is 0. 973832 cnn+eln CA Acuuyracy IS 0. 823473 cnn+eln CA Acuuyracy 2X 0. 981738
cnn+eln. Overall Accuracy is 0. 949029 cnn+eln. Average Accuracy is 0・ 973815 cnn+eln. kappa is 0. 943 1 99
k 打 k WJU-S CS ：3
file 陆畑 l!5*ft Took Oedo^ Wh 如 Help
Testing] me =
20 <0 60 K IK t» »C « 18C 2»
%%Test 1 for PaviaUniversity
%% STEP 1: Initialize
cnnConfig=config ();
[theta, meta] 二cnnInitParams(cnnConfig):
n^class 二{}・dimension;
load PaviaU;
load PaviaU_gt;
[images, labels, indexs, samimage」二loadtrain(paviaU, paviaU_gt);
d 二・ dimension;
images 二reshape (images, d(l), d(2), d(3),[]);
%% STEP 2: Train (The Cnn And Elm)
=1；
=30;
•
f
newtheta 二minFuncSGD(@(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, 900); %% STEP 3: Test [testimages, testLabels, testlndexs]二loadtest(paviaU,
paviaU^gt);
testimages 二 reshape(testImages, d(l),d(2), d(3),[]): [row, col] 二 size(paviaU_gt);
testK 二 cnnExtract(newtheta, testimages, cnnConfig, meta); [TestingTime, TestingAccuracy, testoutputlabel, actualoutputs] 二 elmpredict(testK, testLabels ，, testlndexs, 1, InputWeight, BiasofHid denNeurons,OutputWeight, NumberofOutputNeurons, row, col); predimagel 二zeros(row, col);
predimagel(testindexs)二testoutputlabel; figure, imagesc(predimagel); axis off :
[OA, kappa, AA, CA]= calcError (testLabels, predimagel(testindexs), 1:n_class);
Overall Accuracy is %f\n , OA): Average Accuracy is %f\n , AA); CA Acuuyracy is %f\n‘ , CA);
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 center
cnnConfig 二 config3():
fprintf( cnn+elm fprintf C cnn+elm
fprintf C cnn+elm fprintf C cnn+elm
kappa is %f\n' , kappa);
[theta, meta] 二cnnInitParams(cnnConfig):
n_class 二{}・dimension;
load Pavia;
load Pavia_gt;
[images, labels, indexs] 二loadtrain3(pavia, pavia_gt):
d 二\1;・ dimension;
images 二reshape(images, d(1), d(2), d(3),[]);
%% STEP 2: Train (The Cnn And Elm)
=1；
=25;
•
f
newtheta 二minFuncSGD(@(x, y, z)
cnnCost (x, y, z, cnnConf ig, meta), theta, images, labels, options); K 二cnnExtract(newtheta, images, cnnConfig, meta);
[TrainingTime, TrainingAccuracy, InputWeight, BiasofHiddenNeurons, OutputWeight, NumberofOutputNeurons] =elmtrain(K, labels' , 1, 850); %% STEP 3: Test
[testimages, testLabels, testlndexs] = loadtest3(pavia, pavia_gt); testlmages 二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);
predimagel^zeros(row, col);
predimagel(testlndexs)=testoutputlabel;
figure, imagesc(predimagel);
axis off
[OA, kappa, AA, CA]= calcError (testLabels, predimagel(testlndexs), 1:n_class);
fprintf(' cnn+elm Overall Accuracy is %f\n, OA);
fprintf(' cnn+elm Average Accuracy is %f\n' , AA);
fprintf(' cnn+elm CA Acuuyracy is %f\n‘,CA);
fprin tf ('cnn+elm kappa is %f \n?, kappa);
figure, imagesc(pavia_gt);
axis off
The main code of Salinas of Matlab
%% Convolution Neural Network And Extrem Learning Machine
%%Test 1 for Salinas
%% STEP 1: Initialize
cnnConf ig=conf ig ();
[theta, meta] 二cnnlnitParams(cnnConfig):
n class 二{}・dimension;
load Salinas_corrected;
load Salinas_gt;
[images, labels, indexs, samimage]二loadtrain(salinas_corrected, salinas_gt);
d 二tl;・ dimension;
images 二reshape(images, d(1), d(2), d(3),[]);
%% STEP 2: Train (The Cnn And Elm)
=1；
=30;
9
newtheta 二minFuncSGD(@(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^, labels' , 1, 2300);
%% STEP 3: Test
[testimages, testLabels, testlndexs] 二
loadtest(salinas_corrected, salinas_gt);
testlmages 二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' , testlndexs, 1, InputWeight, BiasofHid denNeurons, OutputWeight, NumberofOutputNeurons, row, col); predimagel二zeros(row, col);
predimagel(testlndexs)二testoutputlabel;
figure, imagesc(predimagel);
axis off:
[OA, kappa, AA, CA]= calcError (testLabels, predimagel(testlndexs),
1:n_class);
fprintf(' cnn+elm Overall Accuracy is %f\n, OA);
fprintf(' cnn+elm Average Accuracy is %f\n' , AA);
fprintf(' cnn+elm CA Acuuyracy is %f\n, CA);
fprintf C cnn+elm kappa is %f\n，, kappa);
figure, imagesc(salinas_gt);
axis off:
六、参考文献
[1]童庆禧，张兵，郑兰芬.高光谱遥感原理.技术与应用[M].北京：高等教育出版社，2006
[2],C・ I. 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. mixture analysis-based compression for hyperspectral image analysis -J]・ IEEE Transanctions on Geoscience and Remote Sensing, 2004, 42 (4): 875-891
[4], of a 3~D hyperspectral instrument for skin imaging
applications [J]・IEEE Transanctions on Instrumentation and Measurement, 2009, 58 (3): 631-638
[5]J・ Li, J. M. Bioucas, A. hyperspectral image segmentation using subspace multinomial logistic regression and markov random fields [J]・IEEE Transanctions on Geoscience and Remote Sensing, 2012, 50(3): 809-823 [6]P. E・ Lehner, L・ Adelman, R・ J. Distasio. Confirmation bias in
the analysis of remote sensing data [J]・ IEEE Transactions on Systems, Man and Cybernetics, 2009, 39(1):218-226・
[7]Bouvrie J・ Notes on Convolutional Neural Networks[J]・ Neural
Nets, 2006.
[8]Ren S, He K, Girshick R, et al. Fas ter R-CNN: Towards Real-Time Object Detection with Region Proposal Networks・[J]・ IEEE
Transactions on Pattern Analysis & Machine Intelligence, 2016:1 T・。