K-means聚类引导的无人机遥感图像阈值分类方法

时间：2024-07-28

白俊龙，王章琼，闫海涛

(1.武汉工程大学土木工程与建筑学院，武汉 430073；2.中交第二公路勘察设计研究院有限公司，武汉 430052)

0 引言

无人机遥感技术以其时效性强、机动灵活、获取成本较低等优势得以快速发展[1]，应用范围也日趋广泛，是传统遥感手段的有力补充[2-3]。近年来，为充分利用无人机遥感影像丰富的信息数据，避免“数据丰富，信息贫乏”[4-5]现象的出现，越来越多的研究者开始了高分辨率遥感图像的处理研究[6]。基于机器学习的图像分类方法是现阶段的研究热点与重点[7]，但此类方法通常需要提供大量的训练样本来保证后续的分类精度[8-9]，且目前各相关领域能够满足模型训练条件的训练样本数据十分稀少，导致该方法还难以用于实际工程应用。而无监督聚类方法不需要训练样本，仅凭地物地磁辐射强弱在遥感图像上所反映的光谱信息，即可将数据按其自然分布特性进行聚类，自动判别地物类别，并能克服人工解译所带来的主观性因素影响[10]。

K-means算法是聚类算法中最为流行的算法，具有思想简单、处理快速及效果好等优点。自1967年MacQueen[11]首次使用后，K-means算法在聚类算法中得到了广泛的应用。但在使用K-means算法时，其聚簇个数k需要事先确定，若k值选取不当会使最终的聚类结果陷入局部最优。而在实际操作中，往往因为数据量过大和缺乏经验导致k值很难确定；且该算法对初始值十分敏感[12]，若k值选取得过小，则会导致同一簇内数据对象差异很大，若k值选择过大，则会导致不同簇间差异很小[13]。故在使用K-means算法进行聚类时，首先应该确定一个最优聚类数目k。

鉴于K-means算法的优点与缺点，结合无人机高分辨率遥感图像的特点，本文提出一种K-means聚类引导的阈值分类方法。首先采用K-means聚类对图像进行聚类处理，实现图像的初始分割，其中利用Average Silhouette指标值来定量确定最优聚类数目k；然后对初始分割结果进行阈值分割和图像优化处理完成对象的提取；再对所有提取结果进行合并实现全图的识别与分类。并基于本文提出的分类方法处理步骤，采用MATLAB/GUI平台，设计了适用于无人机遥感图像分类处理系统。

1 K-means聚类引导的阈值分类方法

1.1 K-means聚类

K-means算法是一种基于划分聚类的无监督学习算法[14]，具有简单、快速、易于解释、性能高，伸缩性好等多项优点[15]。K-means算法通常以欧式距离作为衡量数据对象间相似度的指标，相似度与数据对象间的距离成反比，即对象间距离越小，则相似度越大[16]，其核心思想为[13]：①从给定的数据集中随机选取k个初始聚类中心；②计算其余数据对象与聚类中心的欧氏距离，根据欧式距离，将数据对象分配到与其相距最近的聚类中心所在的簇中；③计算每个簇中数据对象的平均值作为新的聚类中心；④进行下一次迭代，直到聚类中心不再变化或达到最大的迭代次数停止。

1.2 聚类数目k的确定

由此可以看出，聚类个数k值选取是否得当会直接影响最终的聚类效果，故使用K-means算法进行聚类首先需要确定一个合适的k值。传统K-means聚类过程中k值多依靠经验确定，聚类结果容易陷入局部最优。本文利用Average Silhouette指标值来定量确定最优聚类数目k，从而提高初分割精度与处理效率、减少后续手动处理步骤。轮廓图(Silhouette)是一种用来刻画聚类效果的度量，每个点的Silhouette值都表示相比于其他簇，这个点与本簇内的其他点的相似程度。

对于样本点i，其Silhouette值[17]Si的计算公式如下：

(1)

Average Silhouette是计算所有对象的Silhouette值的平均值，反映了分配聚类对象的整体效果，即聚类质量[18]，根据Average Silhouette值可以评估分类结果的准确度或合适度。在聚类数目与Average Silhouette指标值的关系曲线图中，Average Silhouette指标值往往会随聚类数目的增加出现先增大后减小的变化趋势，图中Average Silhouette指标值最大的点所对应的聚类数目，则表示该数据集的最优聚类数目k(图1)。

图1 聚类数目与Average Silhouette值关系示意图Fig.1 Schematic diagram of the relationship between the number of clusters and the Average Silhouette value

1.3 分类步骤

本文提出K-means聚类引导的无人机遥感图像阈值分类方法，主要包括以下步骤：

1)根据已有的无人机遥感图像数据计算Average Silhouette指标值，确定当前图像数据集的最优聚类数目k值；对待处理图像进行K-means聚类初分割；根据K-means聚类结果，目视判定对象类别，将同类别对象进行组合并手工剔除非目标区域，得到新对象。

2)对新对象进行阈值分割，将得到的二值图像进行形态学等优化处理。

3)重复步骤1)和2)直至完成所有区域的处理。

4)最后对处理完的所有区域进行合并，得到完整分类结果。

具体分类流程参见图2。

图2 K-means聚类引导的阈值分类流程Fig.2 The process of threshold classification under K-means clustering guidance

基于MATLAB/GUI平台，对本文提出的分类方法处理步骤进行集成，开发了无人机遥感图像分类处理系统，后续使用该系统完成无人机遥感图像分类的全部处理，系统操作界面如图3所示。

另外，我国有15个崩塌、滑坡和泥石流多发区，分别是：横断山区、黄土高原地区、川北陕南地区、川西北龙门山地区、金沙江中下游地区、川滇交界地区、汉江安康—白河地区、川东大巴山地区、三峡地区、黔西六盘水地区、湘西地区、赣西北地区、赣东北上饶地区、北京怀柔—密云地区、辽东岫岩—凤城地区。

图3 无人机遥感图像处理系统界面Fig.3 UAV remote sensing image processing system interface

2 K-means聚类引导的图像阈值分类

2.1 原始图像与k值确定

采用高精度无人机对桂林至柳州高速公路进行拍摄，将获得的高分辨率可见光遥感图像作为原始数据源，航片的地面分辨率为0.04 m。工作区地势西北高东南低，地表起伏较大，区域内可分为构造侵蚀剥蚀低山丘陵、构造侵蚀剥蚀丘陵、侵蚀堆积河谷阶地、溶蚀堆积孤峰平原4个地貌区。根据《土地利用现状分类》(GB/T 21010—2017)[19]，该区地物类型主要包括林地、草地、耕地、交通运输用地、水域及其他土地等多种土地利用类型。选取该工作区所拍摄的部分图像进行拼接，生成待分类的无人机遥感图像，如图4所示。

图4 拼接之后的无人机遥感图像Fig.4 UAV remote sensing image after stitching

从采用无人机遥感系统获取的高分辨率遥感图像数据来看，针对同一批次的无人机遥感图像数据集，无人机的飞行高度相对固定，单张图像覆盖的面积范围有限，各张遥感图像中所包含的地物类别数量相对较少且基本相等，故处理同一批次的遥感图像数据，只需计算1个k值。通过对本次获取的无人机遥感图像数据集进行大量测试验证，根据式(1)计算对应的Silhouette指标值，并根据聚类数目与Average Silhouette指标值的关系曲线图(图5)，可以确定一个最优聚类数目k为4，故在使用K-means对图像进行聚类时，不会陷入局部最优，提高后续处理效率。

图5 聚类数目与Average Silhouette值的关系Fig.5 Relationship between the number of clusters and the Average Silhouette value

2.2 K-means聚类初分割

首先，使用K-means聚类对原始图像进行初分割。根据2.1节计算确定的最优聚类数目，将初始聚类值设置为4，其初始聚类分割结果如图6所示。通过目视判断可知，初始分割可将部分地物较为完整地分离出来，如图6(b)只包含“林地、草地”(本文将林地和草地按同一类型划分)；绝大部分分割结果包含2类及以上地物，图6(a)包含“其他土地”、“水域”和小部分“交通运输用地”，图6(c)包含“交通运输用地”、“水域”和小部分“其他土地”、图6(d)主要包含“林地、草地”和“水域”。

(a)分割结果1 (b)分割结果2 (c)分割结果3 (d)分割结果4

根据目视判定结果将图6(b)和(d)进行组合，可得到新的处理对象。在系统中选用“抠图”操作将新对象中的“水域”剔除，得到完整的“林地、草地”区域(图7)，完成K-means聚类引导的初分割。

图7 图像组合与抠除Fig.7 Image combination and cutout

2.3 阈值分割与优化

处理得到完整的“林地、草地”区域后，需对图像进行阈值分割和进一步的优化处理，生成二值图像。使用可视化图像阈值选择GUI工具对“林地、草地”进行阈值分割，thresh_tool函数会先自动计算出一个自适应阈值，通常能对图像进行很好地分割；若不满意该结果，也可左右移动灰度直方图中的阈值灰色竖线，调整阈值直至最佳分割效果。如图8所示，对“林地、草地”进行分割的阈值为32。无人机遥感图像分辨率高，在使用K-means聚类对两种相互交融的地物(林地之间的其他土地)进行初分割时，难以得到精准的边界，进行阈值分割之后得到的二值图像往往包含背景噪声。中值滤波是一种非线性操作，使用medfilt2函数对二值图进行中值滤波，可以滤除图中的椒盐噪声，同时保留边缘。再使用imfill函数对图中的密集独立的的孔洞区域进行填充使其连成一整块区域，实现二值图像的优化处理(图9)。

图8 图像阈值选择GUI工具Fig.8 GUI tool for image threshold selection

图9 优化二值图像Fig.9 Optimized binary image

2.4 标记与合并

通过bwlabel函数对二值图像中各个分离部分进行标注，返回一个和二值图像大小相同的标记矩阵L，包含了二值图像中每个连通区域的类别标签。前面经过K-means初分割和阈值分割之后得到的对象中只包含一种地物，故标记矩阵L中只有一种类别标签，其对应的标签值为1。使用label2rgb函数，可将标记矩阵L转换为伪色彩图像，至此即完成“林地、草地”的识别与标记提取。重复上述步骤，分别完成剩下K-means聚类结果的处理，最终从原始图像中识别出“林地、草地”，“水域”，“交通运输用地”和“其他土地”4类地物。

对4类地物标签进行合并，得到彩色标记矩阵，将原始遥感图像与彩色标记矩阵叠加即为最终的图像分类结果。图10为设计的无人机遥感图像处理系统子界面。由图可知，该界面除了对地物标签进行合并与展示最终的分类结果外，还可以计算某类地物的像素面积大小及其占全图总面积的比值。

图10 图像处理系统子界面Fig.10 Sub interface of image processing system

3 结果与分析

3.1 精度评价指标

混淆矩阵分析是进行图像分类精度评价的常用方法[20]。本文使用混淆矩阵计算出的生产者精度、用户精度、总体精度和Kappa系数作为评价指标，对无人机遥感图像分类结果进行精度评价和分析。

1)生产者精度(producer’s accuracy,PA)。表示在分类中某一类别被正确分类的像元数与此类别真实参考像元之间的比率，公式为：

PA=xii/x+i，

(2)

式中：xii为i类地物正确分类的像元数目；x+i为第i列的总像元数量。

2)用户精度(user’s accuracy,UA)。指正确分到某类的像元总数与所有被分为该类的像元总数之间的比率，公式为：

UA=xii/xi+，

(3)

式中，xi+为第i行的总像元数。

3)总体分类精度(overall accuracy,OA)。指被正确分类的像元个数与总像元个数的比值，公式为：

(4)

式中：n为总类别数；N为用于精度评估的像素总数。

4)Kappa系数。用于检验遥感图像分类对于真实地物判断的正确性程度，可以反映所分类别与遥感图像之间的一致性。Kappa系数通常介于0～1之间，Kappa值越大表示分类精度越高。公式为：

(5)

3.2 精度评价

利用449个验证样本点数据，分别建立4种土地利用类型数据集在整体区域的混淆矩阵(表1)，计算出的OA为91.09%，Kappa系数为0.88。可见，使用本文方法对无人机遥感图像分类得到的OA和Kappa系数都很高，说明分类结果与地表真实信息的总体一致性程度高，该方法适用于无人机遥感图像的分类处理。

表1 混淆矩阵Tab.1 Confusion matrix

根据表1可知，用户精度：水域>交通运输用地>林地、草地>其他土地；生产者精度：交通运输用地>水域>林地、草地>其他土地。可以发现其他土地利用类型的UA和PA均为最低，林地、草地次之，分析其原因，该类型分类精度低是因为紧挨林地、草地的其他土地对象中混入部分林地、草地所致；水域和交通运输用地的两种精度都很高，是因为水域多呈面状或带状分布，交通运输用地呈带状分布，区域内地物分布均匀、单一，且这两种形状与其他地物接触通常有较为明显的边界，易于识别区分。