|本文共548字，阅读约需2分钟
1
算法介绍
K-means 是一种常见的无监督学习聚类算法，它通过迭代过程将数据集中的样本划分到不同的簇中。在机器学习和数据挖掘领域，K-means 算法的目标是将 n 个观测值分配到 k 个聚类中，使得每个观测值都属于距离其最近的质心（即集群中心）所在的聚类，同时所有质心是各个聚类中所有点的均值。
在遥感地学领域中，K-means聚类算法被广泛应用于图像分类、土地覆盖分类以及变化检测等任务。该算法通过对遥感图像的像素值进行处理，将每个像素视为一个多维特征向量，然后根据其光谱特征将其分配到预设的类别中。
本文采用Landsat 7影像作为演示数据，利用K-means进行聚类分析，采用CH分数确定最佳聚类数。
2
数据导入
主要基于gdal读取数据，并将其以数组形式进行存储，可视化短波红外波段。
from __future__ import print_function, division# 导入库from osgeo import gdal, gdal_arrayimport numpy as npimport matplotlib.pyplot as plt%matplotlib inlinegdal.UseExceptions() #设置 GDAL 库是否使用异常处理gdal.AllRegister() #注册 GDAL 中的所有驱动程序# 读取影像数据和roiimg_ds = gdal.Open('/home/mw/input/landsat76958/Landsat7/LE70220491999322EDC01_stack.gtif', gdal.GA_ReadOnly)#将数据其存为一个数组img = np.zeros((img_ds.RasterYSize, img_ds.RasterXSize, img_ds.RasterCount),gdal_array.GDALTypeCodeToNumericTypeCode(img_ds.GetRasterBand(1).DataType))for b in range(img.shape[2]):img[:, :, b] = img_ds.GetRasterBand(b + 1).ReadAsArray()# 可视化plt.figure(figsize=(10,8))plt.imshow(img[:, :, 4], cmap=plt.cm.Greys_r)plt.title('SWIR1')
3
确定聚类数
Calinski-Harabasz (CH) 分数，也被称为 Calinski指数或 Variance Ratio Criterion（方差比准则），是一种用于评估聚类效果的内部验证指标。在K-means聚类中，我们经常需要量化不同聚类方案的好坏，以便确定最佳的类别数量。
Calinski-Harabasz分数基于组间和组内的协方差矩阵来计算，通过比较类间离散度与类内离散度的关系，得出一个数值结果。较高的Calinski-Harabasz分数意味着更好的聚类性能，即类间的差异相对于类内的差异更大。
from sklearn.metrics.cluster import calinski_harabasz_scorefrom sklearn.cluster import KMeans# 将所有波段数据堆叠为单个多维特征矩阵，每一行代表一个像素点pixel_vectors = img.reshape((-1, img.shape[2]))# 确保数据是浮点类型，因为KMeans默认处理连续数值pixel_vectors = pixel_vectors.astype(np.float64)#确定聚类数量for i in range(2,7):kmeans=KMeans(n_clusters=i,random_state=123).fit(pixel_vectors)score=calinski_harabasz_score(pixel_vectors,kmeans.labels_)print('聚类%d簇的calinski_harabaz分数为：%f'%(i,score))
4
模型构建
如上所示，聚类4簇的CH最大，因此将其分成4类，并对结果进行可视化展示。
i=4 #聚类4簇的CH最大，因此将其分成4类# 创建并训练KMeans模型estimator = KMeans(n_clusters=i)estimator.fit(pixel_vectors)# 聚类结果是一个标签数组，对应于原始数据中每个像素的类别labels = estimator.labels_# 为了展示分类结果，需要将标签转换回原始图像尺寸image_clustered = labels.reshape(img.shape[:2])cmap = plt.cm.get_cmap('terrain', i)# 可视化聚类结果plt.figure(figsize=(10,8))plt.imshow(image_clustered, cmap=cmap)plt.title('Clustered Image')# 创建图例colors = [cmap(i) for i in range(i)]labels = ['Cluster {}'.format(i) for i in range(i)]plt.legend(handles=[plt.Rectangle((0,0),1,1,fc=color) for color in colors], labels=labels)plt.show()
5
代码获取
公众号后台回复【240313】，即可获取完整的代码及示例数据！
喜欢也行，不喜欢也行；如果觉得有用处的话，还请点点右下角的赞/在看，记得关注我哟！！！
