【机器学习与遥感】sklearn与rasterio实现遥感影像非监督分类

在学习遥感的过程中，我们都了解到了监督分类与非监督分类，二者是遥感解译的基础。之前更多的是使用Erdas与ENVI来进行这两种分类。这里使用python语言，基于机器学习库sklearn与遥感影像处理库rasterio，使用kmeans动态聚类方法实现非监督分类。

非监督分类

非监督分类是以不同影像地物在特征空间中类别特征的差别为依据的一种无先验（已知）类别标准的图像分类，是以集群为理论基础，通过计算机对图像进行集聚统计分析的方法。根据待分类样本特征参数的统计特征，建立决策规则来进行分类。而不需事先知道类别特征。把各样本的空间分布按其相似性分割或合并成一群集，每一群集代表的地物类别，需经实地调查或与已知类型的地物加以比较才能确定。有监督必须有训练集与测试样本。在训练集中找规律，而对测试样本使用这种规律；非监督没有训练集，只有一组数据，在该组数据集内寻找规律。

sklearn

Scikit-learn（以前称为scikits.learn，也称为sklearn）是针对Python 编程语言的免费软件机器学习库 。它具有各种分类，回归和聚类算法，包括支持向量机，随机森林，梯度提升，k均值和DBSCAN，并且旨在与Python数值科学库NumPy和SciPy联合使用。

rasterio

rasterio是一个操作栅格影像的python库，它是在GDAL和Pyproj等库的基础上封装而来， 它的函数通常接受并返回 Numpy ndarrays。 Rasterio 旨在使处理地理空间栅格数据更高效、更有趣。它也可以直接进行影像的读写，包括弹窗渲染显示影像。

K-Means

K-Means 是发现给定数据集的 K 个簇的聚类算法, 之所以称之为 K-均值 是因为它可以发现 K 个不同的簇, 且每个簇的中心采用簇中所含值的均值计算而成。

优点

• 容易理解，聚类效果不错，虽然是局部最优， 但往往局部最优就够了；
• 处理大数据集的时候，该算法可以保证较好的伸缩性；
• 当簇近似高斯分布的时候，效果非常不错；
• 算法复杂度低。

缺点

• K 值需要人为设定，不同 K 值得到的结果不一样；
• 对初始的簇中心敏感，不同选取方式会得到不同结果；
• 对异常值敏感；
• 样本只能归为一类，不适合多分类任务；
• 不适合太离散的分类、样本类别不平衡的分类、非凸形状的分类。

代码实现过程

1. 依赖

import rasterio as rio
from rasterio.plot import show
from sklearn import cluster
import matplotlib.pyplot as plt
import numpy as np
import os
# 处理pyproj安装冲突问题，注意替换为自己的路径
os.environ['PROJ_LIB'] = r'C:\Users\bo\AppData\Local\Programs\Python\Python39\Lib\site-packages\pyproj\proj_dir\share\proj'

需要注意的是os.environ['PROJ_LIB']是用来防止本地pyproj路径冲突的，如果本地无冲突的话，也可以不写。

2.读取原始影像数据

    rs = rio.open(r'E:\测试数据\无人机影像\0926_01.tif')rs_data = rs.read()vmin, vmax = np.nanpercentile(rs_data, (2, 98))

numpy.nanpercentile(arr, q, axis=None, out=None)

• arr : 输入阵列。
• q : 百分数值。
• axis：我们要计算百分位数的轴。否则，它将认为arr是平坦的（在所有轴上工作）。 axis=0意味着沿着列， axis=1意味着沿着行工作。
• out :不同的数组，我们要把结果放在其中。该数组必须具有与预期输出相同的尺寸。

3.数组纬度转换

多波段影像读取出来的数组格式，与聚类分析所使用的格式不一致，需要进行转置。

    # 波段放在了最后一个维度，多维矩阵转置rs_data_trans = rs_data.transpose(1, 2, 0)print(rs_data.shape, rs_data_trans.shape)# 然后直接把矩阵reshape成机器学习中常用的表格形状，将矩阵a变成列数为x，行数不规定的矩阵，具体行数按照总元素个数除列数，均分得到。# 3维转2维rs_data_1d = rs_data_trans.reshape(-1, rs_data_trans.shape[2])print(rs_data_1d.shape)

• 使用 numpy.transpose ()进行变换，其实就是交换了坐标轴，如：x.transpose(1, 2, 0)，其实就是将x第二维度挪到第一维上，第三维移到第二维上，原本的第一维移动到第三维上，最后的shape为：(3，2，2)

• a.reshape(-1, x)则是将矩阵a变成列数为x，行数不规定的矩阵，具体行数按照总元素个数除列数，均分得到。

4.建立模型，进行训练

    # 建立模型# n_clusters：分类簇的数量。cl = cluster.KMeans(n_clusters=5, n_init='auto')  # create an object of the classifierparam = cl.fit(rs_data_1d)  #

KMeans(n_clusters=8, init='k-means++', n_init=10, max_iter=300, tol=0.0001, precompute_distances='auto', verbose=0, random_state=None, copy_x=True, n_jobs=None, algorithm='auto')

• n_clusters：分类簇的数量。

• init：接收待定的string。kmeans++表示该初始化策略选择的初始均值向量之间都距离比较远，它的效果较好；random表示从数据中随机选择K个样本最为初始均值向量；或者提供一个数组，数组的形状为（n_cluster,n_features），该数组作为初始均值向量。

• n_init：用不同的聚类中心初始化值运行算法的次数，最终解是在inertia意义下选出的最优结果；默认值为10。

• max_iter：最大的迭代次数。

• tol：表示算法收敛的阈值。

• precompute_distance：接收Boolean或者auto。表示是否提前计算好样本之间的距离，auto表示如果nsamples*n>12 million，则不提前计算。

• verbose：0表示不输出日志信息；1表示每隔一段时间打印一次日志信息。如果大于1，打印次数频繁。

• random_state：表示随机数生成器的种子。

• n_jobs：表示任务使用CPU数量；若值为 -1，则用所有的CPU进行运算。

5. 提取计算结果

提取计算结果，并将多维数组的转置还原

    # 还原数据img_cl = cl.labels_img_cl = img_cl.reshape(rs_data_trans[:, :, 0].shape)

6. 保存分类结果数据

将监督分类后的影像保存本地

    prof = rs.profile# 更新为1个波段prof.update(count=1)with rio.open(r'E:\测试数据\无人机影像\0926_01result.tif', 'w', **prof) as dst:dst.write(img_cl, 1)

7.结果对别展示

    fig, (ax1, ax2) = plt.subplots(figsize=[15, 15], nrows=1, ncols=2)show(rs, cmap='gray', vmin=vmin, vmax=vmax, ax=ax1)show(img_cl,cmap='gray', ax=ax2)ax1.set_axis_off()ax2.set_axis_off()fig.savefig("pred.png", bbox_inches='tight')plt.show()

对比效果图如下所示：可能我这块准备的影像不太适合

整体流程总结