高光谱图像的SVM分类

#本人小菜鸟一只，目前在学习高光谱图像分类的相关知识，感觉分享的代码实在少得可怜，前几天看了几篇大佬的代码分享，若有所思，按自己看得懂的方式重新整理了一下，分享一下。

#这也是我写的第一篇博客，之后可能也会用这样的方式做分享或者管理我的代码，实在是第一次写，有很多地方不熟悉，如果写的水平太低恶心到你了，实在不好意思。

#在此，也把代码的来源，几位的大佬的博客地址引一下
链接1: xiao韩.链接2: 哈士奇说喵.链接3: 冰机灵.

开始

1.先说一下学习的经过吧，首先是看了链接1中xiao韩老师的关于《基于支持向量机的几种数据预处理的高光谱数据集分类分析》，主要理解原理，就将PCA+LDA+SVM这种组合的代码换用其他的数据跑通了一遍，并未追求高精度的输出结果，仅仅想着把代码跑通，看懂，知道参数怎么修改，可以怎么调试。
本文也主要是对xiao韩老师老师没有说到的地方，加了一些我自己的理解。
2. 对输入图像与结果输出的可视化方面借鉴了哈士奇说喵老师的，主要是利用了spectral的光谱工具包，方便大家能够认识高光谱图像，然后我就发现了冰机灵老师的代码分享，它们的代码风格相似，都各有所长，大家可以综合的看。
3. 本博客要说的内容处理除了是将流程走一遍记录保存代码外，多加了一些自己的理解，和三位老师关于主成分分析（PCA）的降维、线性判别分析（LDA）的降维、网格搜索法（GS）寻找最佳参数在这三个地方没说清的地方做一个小补充，主要以实验为主，不注重原理解讲解。废话不多说上代码！！！

###补充：有需要的同学可以自行下载: 高光谱数据.
我利用的是Salinas的数据，分为两部分Salinas.mat为需要作为输入,训练的图像数据, Salinas_gt.mat为ground_truth,将作为输出的标签数据。

代码

加载遥感图像的.mat数据进行呈图显示

# 加载数据
import matplotlib.pyplot as plt
import numpy as np
from scipy.io import loadmat
import spectral
from functools import reduce # 获取mat格式的数据，loadmat输出的是dict，需要进行定位
input_image = loadmat('E:/Hyperspectral/hyperspectral_code/data/Salinas.mat')['salinas']
output_image = loadmat('E:/Hyperspectral/hyperspectral_code/data/Salinas_gt.mat')['salinas_gt']# input_image.shape      # (512, 217, 224)
# output_image.shape     # (512, 217)
# np.unique(output_image)  # array([ 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16],dtype=uint8)

此处加载数据时需要注意【‘salinas’】表明你要导入的变量的名字，这个地方大小写错了可是会出错的，对于一个新的数据，以防万一，可以使用matlab加载你的数据，然后查看变量的名字，或者修改。

# 统计每类样本所含个数
dict_k = {}
for i in range(output_image.shape[0]):for j in range(output_image.shape[1]):if output_image[i][j] in [m for m in range(1,17)]:#if output_image[i][j] in [1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13]:if output_image[i][j] not in dict_k:dict_k[output_image[i][j]]=0dict_k[output_image[i][j]] +=1print(dict_k)
print(reduce(lambda x,y:x+y,dict_k.values()))# {6: 3959, 7: 3579, 4: 1394, 5: 2678, 15: 7268, 8: 11271, 3: 1976, 2: 3726, 1: 2009, 11: 1068, 12: 1927, 13: 916, 14: 1070, 10: 3278, 9: 6203, 16: 1807}
# 54129

# 展示地物
ground_truth = spectral.imshow(classes = output_image.astype(int),figsize =(9,9))

大佬在这又给出了颜色绘制的变化，可以根据植被，山地等类型绘制适合的颜色，增加图片的可读性。

# 不同颜色绘制
my_color =np.array([[255,255,255],[184,40,99],[74,77,145],[35,102,193],[238,110,105],[117,249,76],[114,251,253],[126,196,59],[234,65,247],[141,79,77],[183,40,99],[0,39,245],[90,196,111],[50,140,100],[70,140,200],[100,150,170]])ground_truth = spectral.imshow(classes = output_image.astype(int),figsize =(9,9),colors=my_color)

重构用到的类，转化为CSV文件

# 除掉 0 这个非分类的类，把所有需要分类的元素提取出来
need_label = np.zeros([output_image.shape[0],output_image.shape[1]])   #(512,217)的全零矩阵
for i in range(output_image.shape[0]):for j in range(output_image.shape[1]):if output_image[i][j] != 0:#if output_image[i][j] in [1,2,3,4,5,6,7,8,9]:need_label[i][j] = output_image[i][j]new_datawithlabel_list = []
for i in range(output_image.shape[0]):for j in range(output_image.shape[1]):if need_label[i][j] != 0:c2l = list(input_image[i][j])c2l.append(need_label[i][j])new_datawithlabel_list.append(c2l)new_datawithlabel_array = np.array(new_datawithlabel_list)
# print(new_datawithlabel_array.shape)  # new_datawithlabel_array.shape (54129,225),包含了数据维度和标签维度，数据225维度，也就是224个波段，最后255列是标签维

from sklearn import preprocessing
data_D = preprocessing.StandardScaler().fit_transform(new_datawithlabel_array[:,:-1])  #(arry[:, :-1])按逆序展示，去掉了最后一列
# print(data_D.shape)  #(54129, 224)
#data_D = preprocessing.MinMaxScaler().fit_transform(new_datawithlabel_array[:,:-1])data_L = new_datawithlabel_array[:,-1]
# print(data_L.shape)  #(54129, 224)# 将结果存档后续处理
import pandas as pd
new = np.column_stack((data_D,data_L))
new_ = pd.DataFrame(new)
new_.to_csv('E:/Hyperspectral/hyperspectral_code/salinas.csv',header=False,index=False)
# 生成csv文件后，就可以直接对该文件进行操作
print('Done')

训练模型并存储模型（PCA-LDA-GS-SVM）

# 导包
import joblib
from sklearn.model_selection import KFold
from sklearn.model_selection import train_test_split
import numpy as np
from sklearn.svm import SVC
from sklearn import metrics
from sklearn import preprocessing
from sklearn.decomposition import RandomizedPCA
import pandas as pd
from sklearn.grid_search import GridSearchCV
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from time import time# 导入数据集切割训练与测试数据data = pd.read_csv('E:/Hyperspectral/hyperspectral_code/salinas.csv',header=None)
data = data.as_matrix()
data_D = data[:,:-1]  #(54129, 224)
data_L = data[:,-1]   #(54129,)
data_train, data_test, label_train, label_test = train_test_split(data_D,data_L,test_size=0.5)# 模型训练与拟合linear  rbf  poly
t0 = time()
pca = RandomizedPCA(n_components = 70, whiten=True).fit(data_D)
X_train_pca = pca.transform(data_train)
X_test_pca = pca.transform(data_test)lda = LinearDiscriminantAnalysis(n_components = 14).fit(X_train_pca,label_train)
X_train_ida = lda.transform(X_train_pca)
X_test_ida = lda.transform(X_test_pca)print('OK')# param_grid = {'C': [1e3, 5e3, 1e4, 5e4, 1e5],
#               'gamma': [0.0001, 0.0005, 0.001, 0.005, 0.01, 0.1], }
# param_grid = {'C': [10, 20, 100, 500, 1e3],
#               'gamma': [0.001, 0.005, 0.01, 0.05, 0.1, 0.125], }
# clf = GridSearchCV(SVC(kernel='linear', class_weight='balanced'), param_grid)clf = SVC(kernel = 'rbf',gamma=0.1,C=20)   # paviau:18-10 rbf 0.1  100  linear 0.125  20  poly 0.1 100  indian:80-40rbf 0.01 100  linear 0.1  20  poly 0.1 100
clf.fit(X_train_ida,label_train)            #salinas:70-14 rbf 0.1  20    linear 0.125  10  poly 0.1 100
# clf.fit(data_train,label_train)
pred = clf.predict(X_test_ida)accuracy = metrics.accuracy_score(label_test, pred)*100
print(accuracy)    # 95.60687234435618
# print(clf.best_estimator_)
print("done in %0.3fs" % (time() - t0))  # done in 13.843s
# 存储结果学习模型，方便之后的调用
joblib.dump(clf, "salinas_MODEL.m")

导入time包的原因是作者在比较几种处理方式时，考虑了处理时间的长短

————————————————————分隔符——————————————————————
分隔符内的内容与主程序段是独立的，主要是我对主成分分析（PCA）的降维、线性判别分析（LDA）的降维、网格搜索法（GS）寻找最佳参数这三个老师没说清的地方进行的补充。
同时，在这里有几个疑问，pca的n_components，lda的n_components怎么确定，原始数据的波段数与降维后的波段数间数据信息量之间含有怎样的关系？作者在这里并没有提到，这是本文在这里会补充的地方。

PCA降维

#经过PCA后，减少计算的复杂程度的基础上，尽可能多的保留原始信息
from sklearn.decomposition import PCA
import numpy as np
from sklearn.preprocessing import StandardScaler
from scipy.io import loadmat
import numpy as np
import matplotlib.pyplot as pltX = loadmat('E:/Hyperspectral/hyperspectral_code/data/Salinas.mat')['salinas']
print (X.shape)
newX = np.reshape(X,(-1,X.shape[2]))
print (newX.shape)
y = loadmat('E:/Hyperspectral/hyperspectral_code/data/Salinas_gt.mat')['salinas_gt']
#pca=PCA( )
pca=PCA(n_components=0.9999)  # n_components可以写数字或小数，此处表示保留99.99%的原始信息
pca.fit(newX,y)
ratio=pca.explained_variance_ratio_
print("pca.components_",pca.components_.shape)
print("pca_var_ratio",pca.explained_variance_ratio_.shape)
#绘制图形plt.plot([i for i in range(X.shape[2])],[np.sum(ratio[:i+1]) for i in range(X.shape[2])])
plt.xticks(np.arange(X.shape[2],step=20))
plt.yticks(np.arange(0,1.01,0.10))
plt.grid()
plt.show()
# 横坐标：表示波段数
# 纵坐标：降维后的所有成分的方差和
# (512, 217, 224)
# (111104, 224)
# pca.components_ (30, 224)
# pca_var_ratio (30,)

此段独立，目的对pca降维的原理代码有所了解，利用了sklearn的工具包，此处pca.components_ (30, 224)表明经PCA变换后，将各个波段按照对主成分信息的贡献成度进行了降序排列，PCA变换后的前30个波段就可以保留原有224个波段的99.99%的主成分信息，其他可能为一些噪声等的无用数据。了解了pca之后，你就自然知道在前一段的n_components应该选择多少波段了 ,文中设置了70个波段，可能是经过测试的结果。

关于LDA的讨论

PCA 和LDA 有很多的相似之处，其本质是将原始的样本映射到维度更低的样本空间中，但是PCA 和LDA 的映射目标不一样： PCA 是为了让映射后的样本具有最大的发散性；而LDA 是为了让映射后的样本有最好的分类性能。所以说， PCA 是一种无监督的降维方法，而LDA 是一种有监督的降维方法。
链接1的文章在处理paviaU的数据时，原有数据的标签类别为9类，n_components 选择了10，而在处理Salinas时，原有数据的标签类别为16类，n_components 选择了10选择了14，猜测是在数据类别的数值附近进行测试的结果。
（话说，我在LDA这个地方的解释是按照我自己的理解，其实写博客本身也就是为了学习，管理、记录代码的，增加对代码的理解程度，万一说错的话希望各位大佬批评指正。）

网格搜索法（GS）寻找最佳参数

# 此段独立，目的利用网格搜索法找到最佳参数gamma和惩罚因子C，将最佳参数可填入上段SVC函数中
# 关于网格搜索法的更多在  https://www.cnblogs.com/wj-1314/p/10422159.html

# 网格搜索法1
from sklearn.datasets import load_iris
from sklearn.svm import SVC
from sklearn.model_selection import train_test_split
import pandas as pd
data = pd.read_csv(r'E:\BaiduNetdiskDownload\hyperspectral_datasets\salinas.csv',header=None)
data = data.as_matrix()
data_D = data[:,:-1]
data_L = data[:,-1]
X_train, X_test, y_train, y_test = train_test_split(data_D,data_L,test_size=0.5)# grid search start
best_score = 0
for gamma in [0.001,0.01,1,10,100]:for c in [0.001,0.01,1,10,100]:# 对于每种参数可能的组合，进行一次训练svm = SVC(gamma=gamma,C=c)svm.fit(X_train,y_train)score = svm.score(X_test,y_test)# 找到表现最好的参数if score > best_score:best_score = scorebest_parameters = {'gamma':gamma,"C":c}print('Best socre:{:.2f}'.format(best_score))
print('Best parameters:{}'.format(best_parameters))

# 网格搜索法2：引入验证集后，为减少偶然性，引入交叉验证
from sklearn.datasets import load_iris
from sklearn.svm import SVC
from sklearn.model_selection import train_test_split,cross_val_score
import pandas as pddata = pd.read_csv(r'E:\BaiduNetdiskDownload\hyperspectral_datasets\salinas.csv',header=None)
data = data.values  #as_matrix
data_D = data[:,:-1]
data_L = data[:,-1]
X_trainval,X_test,y_trainval,y_test = train_test_split(data_D,data_L,test_size=0.2,random_state=0)
X_train ,X_val,y_train,y_val = train_test_split(X_trainval,y_trainval,random_state=1)# grid search start
best_score = 0
for gamma in [0.001]:  #,0.01,1,10,100for c in [1]:     #0.001,0.01,,10,100# 对于每种参数可能的组合，进行一次训练svm = SVC(gamma=gamma,C=c)# 5 折交叉验证scores = cross_val_score(svm,X_trainval,y_trainval,cv=5)score = scores.mean()# 找到表现最好的参数  if score > best_score:best_score = scorebest_parameters = {'gamma':gamma,"C":c}# 使用最佳参数，构建新的模型
svm = SVC(**best_parameters)# 使用训练集和验证集进行训练 more data always resultd in good performance
svm.fit(X_trainval,y_trainval)# evalyation 模型评估
test_score = svm.score(X_test,y_test)print('Best socre:{:.2f}'.format(best_score))
print('Best parameters:{}'.format(best_parameters))
print('Best score on test set:{:.2f}'.format(test_score))

————————————————————分隔符——————————————————————
通过前面我们经过训练后得到了所构造的模型，后面将调用此模型对数据进行预测。后面的这部分内容是学习链接二中哈士奇说喵老师的代码进行了部分的改动。改动的部分主要是因为模型是经过了PCA与LDA的，所以我们测试时也要经过相同的处理步骤。

模型预测在图中标记类

import matplotlib.pyplot as plt
import numpy as np
from scipy.io import loadmat
import spectral# mat文件的导入
import matplotlib.pyplot as plt
import numpy as np
from scipy.io import loadmat
import spectralinput_image = loadmat('E:/Hyperspectral/hyperspectral_code/data/Salinas.mat')['salinas']
output_image = loadmat('E:/Hyperspectral/hyperspectral_code/data/Salinas_gt.mat')['salinas_gt']testdata = np.genfromtxt('E:/Hyperspectral/hyperspectral_code/salinas.csv',delimiter=',')
data_test = testdata[:,:-1]
label_test = testdata[:,-1]#在利用训练的model.m模型进行预测时，因为模型经过了PCA与LDA,此时测试时也要经过相同的处理步骤
pca = RandomizedPCA(n_components = 70, whiten=True).fit(data_test)
data_test_pca = pca.transform(data_test)lda = LinearDiscriminantAnalysis(n_components = 14).fit(data_test_pca,label_test)
data_test_pca_lda = lda.transform(data_test_pca)
# print(data_test_pca_lda.shape)print('ok')# salinas_MODEL.m
clf = joblib.load("salinas_MODEL.m")predict_label = clf.predict(data_test_pca_lda)accuracy = metrics.accuracy_score(label_test, predict_label)*100print(accuracy) # 79.45833102403518# 将预测的结果匹配到图像中
new_show = np.zeros((output_image.shape[0],output_image.shape[1]))
k = 0
for i in range(output_image.shape[0]):for j in range(output_image.shape[1]):if output_image[i][j] != 0 :new_show[i][j] = predict_label[k]k +=1 # print new_show.shape

# 将预测的结果匹配到图像中
new_show = np.zeros((output_image.shape[0],output_image.shape[1]))
k = 0
for i in range(output_image.shape[0]):for j in range(output_image.shape[1]):if output_image[i][j] != 0 :new_show[i][j] = predict_label[k]k +=1 # print new_show.shape# 展示地物
ground_truth = spectral.imshow(classes = output_image.astype(int),figsize =(9,9))
ground_predict = spectral.imshow(classes = new_show.astype(int), figsize =(9,9))

在这里，模型训练测试时的精确度约为95.61%，调用模型是精确度约为79.46%，从预测的图像中也能看出预测的图像中存在错分的现象，但整体来时RSV已经是一个比较良好的分类模型了，经过PCA与LDA的两次降维变换主要是为了在保证准确率的基础上，减少计算的复杂程度，这个思维给我们提供了一个创新的思路。
文中如果有说的不对的地方，请批评指正吧，大家一起来学习吧。