《机器学习》课程项目报告高光谱图像分类——基于CNN和ELM学院信息工程学院专业电子与通信工程学号 35学生姓名曹发贤同组学生陈惠明、陈涛硕士导师杨志景2016 年 11 月一、项目意义与价值高光谱遥感技术起源于 20 世纪 80年代初，是在多光谱遥感技术基础之上发展起来的[1]。
高光谱遥感能够通过成像光谱仪在可见光、近红外、短波红外、中红外等电磁波谱范围获取近似连续的光谱曲线，将表征地物几何位置关系的空间信息与表征地物属性特征的光谱信息有机地融合在了一起，使得提取地物的细节信息成为可能。
随着新型成像光谱仪的光谱分辨率的提高，人们对相关地物的光谱属性特征的了解也不断深入，许多隐藏在狭窄光谱范围内的地物特性逐渐被人们所发现，这些因素大大加速了遥感技术的发展，使高光谱遥感成为 21 世纪遥感技术领域重要的研究方向之一。
在将高光谱数据应用于各领域之前，必须进行必要的数据处理。
常用的数据处理技术方法包括：数据降维、目标检测、变化检测等。
其中，分类是遥感数据处理中比较重要的环节，分类结果不但直接提取了影像数据有效信息，可以直接运用于实际需求中，同时也是实现各种应用的前提，为后续应用提供有用的数据信息和技术支持，如为目标检测提供先验信息、为解混合提供端元信息等。
相对于多光谱遥感而言，由于高光谱遥感的波谱覆盖范围较宽，因此我们可以根据需要选择特定的波段来突显地物特征，从而能够精确地处理地物的光谱信[2]。
目前，许多国家开展大量的科研项目对高光谱遥感进行研究，研制出许多不同类型的成像光谱仪。
高光谱遥感正逐步从地面遥感发展到航空遥感和航天遥感，并在地图绘制、资源勘探、农作物监测、精细农业、海洋环境监测等领域发挥重要的作用。
高光谱遥感技术虽然是遥感领域的新技术，但是高光谱图像的分类一直制约着高光谱遥感的应用[3,4]，因此对其进行研究显得尤为重要。
高光谱遥感图像较高的光谱分辨率给传统的图像分类识别算法提出严峻的挑战。
波段维数的增加不仅加重了数据的存储与传输的负担，同时也加剧了数据处理过程的复杂性，并且由于波段与波段间存在着大量的冗余信息，从而使得传统图像分类算法并不适用于高光谱遥感图像的分类。
传统的分类方法往往需要很多数目的已知类别的训练样本，从而导致计算量大，时间复杂度高。
另外，如果训练样本的数目较少，采用传统分类算法进行分类时分类精度往往是很低的，因此提高分类精度并减少运算量成为高光谱领域的热点问题。
高光谱遥感图像的波段数目多，并且波段与波段间存在着很大的相关性，因此在进行遥感图像的处理（例如分类）之前通常需要进行降维预处理，这样做不仅可以降低数据空间的维数，减少冗余信息，而且还有利于人工图像解译和后续分类处理和地物识别，从而为解决高光谱遥感分类的难点提供了方便[5]。
二、高光谱图像分类的发展与现状高光谱图像分类作为高光谱图像的基础研究，一直是高光谱图像重要的信息获取手段，它的主要目标是根据待测地物的空间几何信息与光谱信息将图像中的每个像素划分为不同的类别。
高光谱图像分类按照是否有已知类别的训练样本的参与，高光谱图像的分类方式分为监督分类与非监督分类[6]。
在遥感图像自动分类中，传统的基于数理统计的分类方法，主要包括最小距离分类、最大似然分类、波谱角分类、混合距离法分类等，主要依赖地物的光谱属性，基于单个像元进行分类。
统计模式识别方法本身的不足：1、最大似然法计算强度大，且要求数据服从正态分布2、K-means聚类分类精度低，分类精度依赖于初始聚类中心3、最小距离法没有考虑各类别的协方差矩阵，对训练样本数目要求低近年来对于神经网络分类方法的研究相当活跃。
它区别于传统的分类方法：在处理模式分类问题时，并不基于某个假定的概率分布，在无监督分类中，从特征空间到模式空间的映射是通过网络自组织完成的。
在监督分类中，网络通过对训练样本的学习，获得权值，形成分类器，且具备容错性。
人工神经网络（ANN）分类方法一般可以获得更高精度的分类结果，因此ANN方法在遥感分类中被广泛应用，特别是对于复杂类型的地物类型分类，ANN方法显示了其优越性。
专家系统分类法也在遥感分类取得了一定的应用。
专家系统是模拟人类逻辑思维的智能系统，将其应用于遥感分类最大的优点就是可以充分利用更多的辅助分类数据。
不过由于专家系统知识库难以建立，影响了它的进一步发展。
支持向量机（SVM）具有严格的理论基础，能较好地解决小样本、非线性、高维数等问题，被成功地应用到多光谱、高光谱遥感图像分类领域。
对于高光谱数据而言，由于波段多、数据量大、数据不确定性等，易受Hughes现象（即训练样本固定时，分类精度随特征维数的增加而下降）影响。
而样本的获取在高光谱分类中往往是一项比较困难的工作，特别是采用高维特征向量时要求每类的样本数都要比特征维数高，因此在高维信息处理中的精度与效率和高光谱遥感信息精细光谱与大数据量之间仍然存在着极大的矛盾。
三、卷积神经网络理论基础卷积神经网络是人工神经网络的一种，它的权值共享网络结构使之更类似于生物神经网络，降低了网络模型的复杂度，减少的权值的数量以节约训练和测试的计算时间。
该优点在网络的输入是多维图像时表现得更加明显，使图像可以直接作为网络的输入，避免了传统识别算法中复杂的特征提取和数据重建过程。
卷积神经网络是为识别二维数据而专门设计的一个多层感知机，其网络对平移、比例变化和倾斜等具有高度不变性[7]。
在CNN中，图像的一小部分（局部感受区域）作为层级结构的最低层的输入，信息再依次传输到不同的层，每层通过一个数字滤波器去获得观测数据的最显着的特征。
这个方法能够获取对平移、缩放和旋转不变的观测数据的显着特征，因为图像的局部感受区域允许神经元或者处理单元可以访问到最基础的特征，例如定向边缘或者角点。
卷积神经网络是一个多层的神经网络，每层由多个二维平面组成，每个平面又由多个独立的神经元组成。
图2为卷积神经网络的整体结构图。
一般地，C层（卷积层）为特征提取层，每个神经元的输入与前一层的局部感受野相连，并提取该局部的特征，一旦该局部特征被提取后，它与其他特征间的位置关系也随之确定下来；S层（下采样层）是特征映射层，网络的每个计算层由多个特征映射组成，每个特征映射为一个平面，平面上所有神经元的权值相等。
特征映射结构采用sigmoid函数等映射函数作为卷积网络的激活函数，使得特征映射具有位移不变性。
此外，由于一个映射面上的神经元共享权值，因而减少了网络自由参数的个数，降低了网络参数选择的复杂度。
卷积神经网络中的每一个特征提取层（C-层）都紧跟着一个用来求局部平均与二次提取的计算层（S-层），这种特有的两次特征提取结构使网络在识别时对输入样本有较高的畸变容忍能力。
卷积神经网络采用有监督学习的方式进行训练，即任何一个训练样本的类别是已知的，训练样本在空间中的分布不再是依据其自然分布倾向来划分，而是根据同类样本和不同类样本中的空间关系进行划分，这需要不断调整网络模型的参数用以划分样本空间的分类边界的位置，是一个耗时且复杂的学习训练过程[8]。
神经网络在进行训练时，所有的网络权值都用一些不同的小随机数进行初始化，这些小的随机数能偶保证网络不会因为权值过大而进入饱和状态，导致训练失败。
神经网络训练算法包括4个主要部分：（1）样本集中取出样本(X,yp )并将其输入网络，X代表图像数组，yp代表其类别；（2）计算此次输入相应的实际输出O p ,这是一个前向传播的过程；（3）用一个指定的损失函数计算出实际输出O p与理想输出Y p的误差；（4）按极小化误差的方法反向传播调整网络权值。
四、极限学习机极限学习机（extreme learning machine）ELM是一种简单易用、有效的单隐层前馈神经网络SLFNs学习算法。
2004年由南洋理工大学黄广斌副教授提出。
传统的神经网络学习算法（如BP算法）需要人为设置大量的网络训练参数，并且很容易产生局部最优解。
极限学习机只需要设置网络的隐层节点个数，在算法执行过程中不需要调整网络的输入权值以及隐元的偏置，并且产生唯一的最优解，因此具有学习速度快且泛化性能好的优点。
极限学习机的网络训练模型采用前向单隐层结构。
设,,m M n分别为网络输入层、隐含层和输出层的节点数，()g x是隐层神经元的激活函数，i b 为阈值。
设有N 个不同样本(),i i x t ，1i N ≤≤ ，其中[][]1212,,...,,,,...,T Tm n i i i im i i i in x x x x R t t t t R =∈=∈ ，则极限学习机的网络训练模型如图1所示。
图1 极限学习机的网络训练模型极限学习机的网络模型可用数学表达式表示如下：()1,1,2,...,Mi ii i j i g x b o j N βω=+==∑式中，[]12,,...,i i i mi ωωωω= 表示连接网络输入层节点与第i 个隐层节点的输入权值向量；[]12,,...,Ti i i in ββββ= 表示连接第i 个隐层节点与网络输出层节点的输出权值向量；[]12,,...,Ti i i in o o o o = 表示网络输出值。
极限学习机的代价函数E 可表示为()1,Nj j j E S o t β==-∑式中，(),,1,2,...,i i s b i M ω== ，包含了网络输入权值及隐层节点阈值。
Huang 等指出极限学习机的悬链目标就是寻求最优的S ，β，使得网络输出值与对应实际值误差最小，即()()min ,E S β 。
()()min ,E S β可进一步写为()()()111,,min ,min ,...,,,...,,,...,i i M M N b E S H b b x x T ωββωωβ=-式中，H 表示网络关于样本的隐层输出矩阵，β表示输出权值矩阵，T 表示样本集的目标值矩阵，H ，β，T 分别定义如下：()()()()()111111111,...,,,...,,,...,M M M M N N m N M N M g x b g x b H b b x x g x b g x b ωωωωωω⨯++⎡⎤⎢⎥=⎢⎥⎢⎥++⎣⎦ 11,T T T T M N M N N Nt T t βββ⨯⨯⎡⎤⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦极限学习机的网络训练过程可归结为一个非线性优化问题。
当网络隐层节点的激活函数无限可微时，网络的输入权值和隐层节点阈值可随机赋值，此时矩阵H 为一常数矩阵，极限学习机的学习过程可等价为求取线性系统H T β= 最小范数的最小二乘解ˆβ ，其计算式为ˆH T β+= 式中H +时矩阵H 的MP 广义逆。
五、 具体实现及主要代码1、训练的样本及其样本图2、原图及进行分类后的图和各个高光谱数据集每一类的分类精度The result of Pavia UniversityThe result of Pavia CenterThe result of Salinas3、主要代码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 = {}.dimension;load PaviaU;load PaviaU_gt;[images, labels, indexs,samimage] = loadtrain(paviaU, paviaU_gt);d = {1}.dimension;images = reshape(images,d(1),d(2),d(3),[]);%% STEP 2: Train (The Cnn And Elm)= 1;= 30;= ;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, testIndexs] = loadtest(paviaU,paviaU_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,BiasofHid denNeurons,OutputWeight,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 %f\n',OA);fprintf('cnn+elm Average Accuracy is %f\n',AA);fprintf('cnn+elm CA Acuuyracy is %f\n',CA);fprintf('cnn+elm kappa is %f\n',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 = {}.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;= ;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, 850); %% STEP 3: Test[testImages, testLabels, testIndexs] = loadtest3(pavia, pavia_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=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 %f\n',OA);fprintf('cnn+elm Average Accuracy is %f\n',AA);fprintf('cnn+elm CA Acuuyracy is %f\n',CA);fprintf('cnn+elm kappa is %f\n',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 = {}.dimension;load Salinas_corrected;load Salinas_gt;[images, labels, indexs,samimage] = loadtrain(salinas_corrected, salinas_gt);d = {1}.dimension;images = reshape(images,d(1),d(2),d(3),[]);%% STEP 2: Train (The Cnn And Elm)= 1;= 30;= ;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, 2300);%% STEP 3: Test[testImages, testLabels, testIndexs] =loadtest(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,BiasofHid denNeurons,OutputWeight,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 %f\n',OA);fprintf('cnn+elm Average Accuracy is %f\n',AA);fprintf('cnn+elm CA Acuuyracy is %f\n',CA);fprintf('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. Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks.[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2016:1-1.。
