深度学习+地理信息小白笔记（matlab实现地物分类）

1 对机器学习的个人理解

机器学习是当今科研领域最为热门的方向之一。它是人工智能的一个分支，简单来说它就是通过算法使得机器能从大量数据中学习规律，从而对新的样本做出智能识别或对未来做出预测。机器学习的核心就是“大量数据”——“特征”——“分类或预测”。

而深度学习则是机器学习的一个新的方向。它的特征在于模拟人体神经元的作用机制来对算法进行优化。近年来它在图像识别与检索、语言信息处理、语音识别等众多领域都取得较为成功的发展。尤其是在计算机视觉领域深度学习一直以来都是重要的研究工具。我们日常所说的“神经网络”，正是深度学习的主要内容。

针对地理信息科学领域而言，机器学习等应用较为广泛。例如土地利用分类这一分类问题，恰恰符合了机器学习的实现过程。土地遥感影像的每个栅格可以视作一个样本，每个栅格的不同波段值就是样本在不同维度的特征值；由于一幅遥感图像由很多像素点构成，这样的多维度矩阵就为机器学习提供了大量的数据以建立模型。

在学习和实践过程中，我们也了解到了一些图像分类的方法。在ENVI中软件提供了KNN算法、支持向量机（SVM）算法和最大似然（朴素贝叶斯）算法等；上述算法都是监督学习方法，即模型需要训练集给定特征以及对应的标签，然后才能对未知标签的样本进行计算。而常用的非监督学习方法有K-means算法，层次聚类算法，主成分分析等，它是根据数据的特征直接将其分为几类而不需要事先学习。而深度学习是在监督学习方法中的神经网络（Neural Network，NN）上发展而来的，具体就是通过人工神经网络（Artificial Neural Network，ANN）为架构，输入各类信号，直接输出期望结果的一种方式，当然深度学习有属于监督学习的方法（如卷积神经网络CNN、多层感知机MLP），也有属于非监督学习的方法（如生成对抗网络GAN）等等。当然除此之外还有一些其他的算法也属于机器学习，如在地理学研究中常见的时间序列分析，相关性分析（空间自相关）等等都是定量化处理地理数据的重要手段。在这些算法的基础上，也衍生出显著性检验、精度验证、数据集创建等一系列相关问题。下文将简单介绍我在学习过程中的一些实践经验。

2 遥感地物分类的代码尝试实现

ENVI中主要提供了最小距离（KNN）、支持向量机（SVM）、最大似然法（朴素贝叶斯）等方法。拟使用Matlab实现上述机器学习的分类算法。

在机器学习中十分重要的是数据集的构建。如果数据集不能正确反映数据整体的分布和分类特征，那么机器学习训练后的模型显然不能得出正确结论。在数据集的构建中，训练集、验证集和测试集的划分十分重要。训练集自然是用于模型的训练，而验证集可以用来评估模型性能，即“调参”，寻找出更合适的模型配置，实际上验证集起到了反馈的作用来调整使得模型在训练集上最优。因此，验证集与训练集往往来源于同一数据序列，即“已知的数据”，由此引出了创建验证集的方法，如k折交叉验证等。在模型训练完成之后就可以使用测试集数据对模型进行最终评估。测试集的数据与训练集或验证集数据不同，属于全新数据，这样可以客观检测模型的实际应用价值。

上述数据集的划分可能存在的问题是欠拟合与过拟合。所谓欠拟合是指模型在训练集上的误差仍然很高，即模型还没有训练完成，而过拟合是指模型在训练集上误差和测试集上误差之间的差距太大，即模型在训练集上表现很好，但在测试集上却表现很差。模型没有理解数据背后的规律，“泛化能力”差——不仅希望模型在训练集上误差小，还希望它在测试集上有良好表现。

图中，训练误差指模型在训练集上的误差随着模型的训练而减小；泛化误差则指模型在测试集上的误差会随着训练先减小后增大，原因正在于一味地重复训练集和提升验证集比例，虽然会使得模型更好地反馈训练集特征，但无法突出整组数据的特征，使得在全新的测试集上表现不良。对此，在数据来源无法获得更多的情况下（遥感影像数据获取与打标签需要耗费大量人工），可以采用以下几种方式来尽量降低过拟合的发生。

①k折交叉验证，前文提到了构建验证集的方法，通过把训练集拆分成k个子数据集来交叉验证提升模型的泛化能力；

②主成分分析对数据进行特征提取降维，减少样本的特征数；

③正则化处理，包括L1（Lasso）和L2（Ridge），简单理解就是在原有损失函数计算的基础上增加参数本身，增加一个惩罚函数来对可能出现的系数过大进行校正；

④集成学习（Bagging/Boosting）也是解决模型过拟合的方案之一。Bagging指的是Bootstrap aggregating，即“有放回抽样”+“聚合”：先对样本进行有放回抽样，对每次抽样接过进行模型训练，最终的测试结果取这些小模型的平均/投票等；Boosting则是在训练过程中不断调整样本权重值，将多个弱分类器组装成一个强分类器。它和Bagging的区别在于每一轮的训练集不变，只是训练集中每个样例在分类器中的权重发生变化。上一轮出错则权重增加。基于上述两种训练模型方式，研究者们先后提出了随机森林（Bagging+决策树）、XGBoost（Boosting+决策树）等方案。

对于遥感地物分类方向，采用监督分类算法训练集需要人工选取ROI并打标签来创建。对此我采用了ENVI中的ROI工具进行创建。而测试集一般来源于分辨率更高的图像或者人工目视解译等.

原始数据是上海市的Landsat8卫星7波段遥感影像，空间分辨率30米，共1001*1001个像元。

在ROI选取中共选择了植被、不透水面、水体三种典型地物采样，共选取16092个样本。

导入Matlab进行数据前期处理。

clear;
clc;
sample2013=csvread("land2013.csv");
sample2017=csvread("land2017.csv");
d2013=imread("20130829.tif");
[m,n,b]=size(d2013);
data2013=reshape(d2013,[m*n,b]);
d2017=imread("20170402.tif");
data2017=reshape(d2017,[m*n,b]);
data2013=double(data2013);
xtrain2013=sample2013(:,3:9);
ytrain2013=sample2013(:,10);
xtrain2017=sample2017(:,3:9);
ytrain2017=sample2017(:,10);
writematrix(xtrain2013,'xtrain2013.csv');
writematrix(ytrain2013,'ytrain2013.csv');
writematrix(data2013,'data2013.csv');

将数据保存后分别用KNN、SVM、Bayes三种方法进行地物分类。

KNN算法没有采用pdist函数计算，而是使用了knnsearch函数。

Mdl=KDTreeSearcher(xtrain2013);
[xn1,~]=knnsearch(Mdl,double(data2013),'k',3);
%确定K个最近邻居中的大多数类别
classnum1=mode(xn1,2);
for i = classnum1
    xi=sample2013(i,10);
end
class_2_mat1=reshape(xi,[m,n]);
imagesc(class_2_mat1)
saveas(1,'knn.jpg')

结果如图所示。

对于SVM算法实现起来有一定难度，原因在于SVM主要用于二分类问题，对于多分类问题需要考虑怎样构造向量机。一般来说有：

①一对多法(OVR -SVM），即训练时把某个类别的样本归为一类，其他剩余的样本归为另一类。这样k个类别的样本就构造出了k个SVM，分类时将未知样本分类为具有最大分类函数值的那类。

②一对一法（OVO-SVM），即在任意两类样本之间设计一个SVM，k个类别的样本就需要设计k*（k-1）/2个SVM。对一个未知样本进行分类时，最后得票最多的类别即为该未知样本的类别。

在Matlab中所提供的函数fitsvm主要用于二分类，不太适合遥感图像的地物分类；可以采用libsvm函数进行分析。过程略有复杂，所以同时采用了python的方法来进行分析。

关于libsvm包，可以参考CSDN中另外博客。

%SVM分类
cmd=["-v",3];
model2=libsvmtrain(ytrain2013,xtrain2013,cmd);
plabel=ones(1002001,1);
data2013=double(data2013);
[predict_label_train,accuracy_train,~]=libsvmpredict(plabel,data2013,model2);
class_2_mat2=reshape(predict_label_train,[m,n]);
imagesc(class_2_mat2)
y_predict=my_MultiSvm(xtrain2013,ytrain2013,data2013);
data=csvread("mat.csv",1,1);
imagesc(data.')
saveas(1,'svm.jpg')

可以看到SVM的分类似更精确，一些在KNN算法中出现的像素点的噪声得到了抑制。

接下来是最大似然法（bayes）的分类方法。最大似然法主要是依据朴素概率原理来对每个像素点是某一类的概率进行预测。

%从训练集的标签数据中找到每个类别的位置信息，以及所占的百分比
%tag 第一列为label中类别的编号，第二列为该类别的像素的数量，第三列为所占的百分比
tag=tabulate(ytrain2013(:));
%将每一个类别的像素的特征值保存到一个类别的矩阵中
x1=xtrain2013(:,1);x2=xtrain2013(:,2);x3=xtrain2013(:,3);x4=xtrain2013(:,4);
x5=xtrain2013(:,5);x6=xtrain2013(:,6);x7=xtrain2013(:,7);
V1(:,1)=x1(ytrain2013==tag(1,1));V1(:,2)=x2(ytrain2013==tag(1,1));
V1(:,3)=x3(ytrain2013==tag(1,1));V1(:,4)=x4(ytrain2013==tag(1,1));
V1(:,5)=x5(ytrain2013==tag(1,1));V1(:,6)=x6(ytrain2013==tag(1,1));
V1(:,7)=x7(ytrain2013==tag(1,1));V2(:,1)=x1(ytrain2013==tag(2,1));
V2(:,2)=x2(ytrain2013==tag(2,1));V2(:,3)=x3(ytrain2013==tag(2,1));
V2(:,4)=x4(ytrain2013==tag(2,1));V2(:,5)=x5(ytrain2013==tag(2,1));
V2(:,6)=x6(ytrain2013==tag(2,1));V2(:,7)=x7(ytrain2013==tag(2,1));
V3(:,1)=x1(ytrain2013==tag(3,1));V3(:,2)=x2(ytrain2013==tag(3,1));
V3(:,3)=x3(ytrain2013==tag(3,1));V3(:,4)=x4(ytrain2013==tag(3,1));
V3(:,5)=x5(ytrain2013==tag(3,1));V3(:,6)=x6(ytrain2013==tag(3,1));
V3(:,7)=x7(ytrain2013==tag(3,1));
%计算判别函数
%计算均值,3×1的列向量
u1=mean(V1)';u2 = mean(V2)';u3 = mean(V3)';
%计算协方差矩阵
c1=cov(V1);c2=cov(V2);c3=cov(V3);
%假设协方差阵相同
c=c1+c2+c3;
%各个类别的判别函数
tic
%先验概率不同
f1=Jug_Func(u1,c,tag(1,3));
f2=Jug_Func(u2,c,tag(2,3));
f3=Jug_Func(u3,c,tag(3,3));
fprintf('生成判别函数用时：%f 秒\n',toc);
%% 把每一个像素都代入判别式中进行计算
%用来保存分类结果
Res=zeros(1001,1001);
tic
rs=zeros(1002001,3);
rs(:,1)=data2013(:,1)*f1(1)+data2013(:,2)*f1(2)+data2013(:,3)*f1(3)+data2013(:,4)*f1(4) ...
+data2013(:,5)*f1(5)+data2013(:,6)*f1(6)+data2013(:,7)*f1(7)+ones(1002001,1)*f1(8);
rs(:,2)=data2013(:,1)*f2(1)+data2013(:,2)*f2(2)+data2013(:,3)*f2(3)+data2013(:,4)*f2(4) ...
+data2013(:,5)*f2(5)+data2013(:,6)*f2(6)+data2013(:,7)*f2(7)+ones(1002001,1)*f2(8);
rs(:,3)=data2013(:,1)*f3(1)+data2013(:,2)*f3(2)+data2013(:,3)*f3(3)+data2013(:,4)*f3(4) ...
+data2013(:,5)*f3(5)+data2013(:,6)*f3(6)+data2013(:,7)*f3(7)+ones(1002001,1)*f3(8);
[~,index]=max(rs,[],2);
for i = 1:3
   Res(index==i)=tag(i,1); 
end
fprintf('分类过程用时： %f 秒\n',toc);
imagesc(Res);
saveas(1,'bayes.jpg')
function result=Jug_Func(u,C,Pw)
    wi=inv(C)*u;
    wio=-0.5*u'*inv(C)*u-0.5*log(det(C))+log(Pw);
    result=[wi(1),wi(2),wi(3),wi(4),wi(5),wi(6),wi(7),wio];
end

最大似然的结果相对来说较为粗糙，原因就在于是概率函数得到的结果，受到人为（选ROI区域）或客观（遥感图像）影响较大，同时也需要训练集尽量大并满足正态分布。

总结来说基本实现了这些功能的代码实现，但仍较为粗糙。事实上了解到python在机器学习特别是深度学习领域运用较为深刻，在之后可以尝试实现这些功能的python代码。

3 机器学习的精度验证

在完成模型的训练后要进行精度验证，特别是遥感地物分类识别需要计算各类精度指标等。常见的有准确度、精确率、召回率、F1、混淆矩阵、Kappa系数、总体精度OA、制图精度PA、用户精度UA、FoM指数等。上述指标基本上都是基于假正例FP（把假的分到真的）、假负例FN（把假的分到别的）、真正例TP（分对了）、真负例TN（把真的分错到别的）这四种在分类中可能出现的情况来计算的。

准确度ACC：(TP+FN)/(TP+FP+FN+TP)

精确率Precision，指识别出来的目标有多少是识别准确的，衡量的是识别结果的正确率：TP/(TP+FP)；

召回率Recall，指所有正确的目标有多少被识别出来了，衡量的是识别结果的查全率：TP/(TP+FN)；

F1，指精确率和召回率的平均：2*precision*recall/（precision + recall）；

总体精度OA：(TP+TN)/(TP+FP+FN+TP)；

制图精度PA：TP/(TP+TN)；

用户精度OA：TP/(TP+FP)；

Kappa系数：基于混淆矩阵的指标，取值为-1到1之间，减小由于样本种类比例不均衡所导致的ACC值偏高的情况。

FoM指数：在模拟未来土地利用变换中FoM指数也是重要指标之一。它考虑了未来土地利用变换情况，增强了对土地利用模拟的模型评价全面性。

其中，A是实际状态为变化而模拟结果为不变的样本个数；B为实际状态为变化且模拟结果为变化的样本个数；C为实际状态为变化且模拟结果为变化，但变化类型有误的样本个数；D为实际状态为不变而模拟结果为变化的样本个数。

尝试代码实现了混淆矩阵相关指标代码与FoM指数代码如下。

%% 混淆矩阵相关函数
function [OA,UA,PA,kappa,precision,recall]=f1(label,predict)
   M=confusionmat(label,predict);
   precision=diag(M)./(sum(M,2));
   recall=diag(M)./(sum(M,1)');
   UA=mean(precision);
   PA=mean(recall);
   OA=sum(diag(M))./(sum(M(:)));
   po=trace(M)/sum(M(:));
   pe=sum((sum(M,2)).*(sum(M,1)'))/(sum(M(:)))^2;
   kappa=(po-pe)/(1-pe);
end
%% FOM计算函数
function [fom,A,B,C,D]=f2(source,real,simu)
    [A,B,C,D]=deal(0);
    siz1=size(source);
    len=siz1(1,1);
    for i=1:len
        if source(i)~=real(i) && simu(i)==source(i) && simu(i)~=real(i)
            A=A+1;
        elseif source(i)~=real(i) && simu(i)~=source(i) && simu(i)==real(i)
            B=B+1;
        elseif source(i)~=real(i) && simu(i)~=source(i) && simu(i)~=real(i)
            C=C+1;
        elseif source(i)==real(i) && simu(i)~=source(i) && simu(i)~=real(i)
            D=D+1;
        end
    end
    fom=B/(A+B+C+D);
end