基于谱聚类算法的三维激光点云数据分类研究

吴翔，王凤艳，林楠，王明常

1.吉林大学 地球探测科学与技术学院，长春 130026；2.吉林建筑大学 测绘与勘查工程学院，长春 130018

0 引言

三维激光扫描技术是通过发射激光来获取物体表面的三维空间信息和物理信息，实现了测绘领域的又一次技术革新[1]。与一些传统的全站仪测量技术、GPS测量技术相比，三维激光扫描技术有着快速、准确和非接触性等特点，广泛应用于城市建设、工程设计和地貌分析等领域[2]。但三维激光扫描技术获取的点云数据量大且杂乱无序，不利于后续的三维场景识别和地物分类等分析与应用。因此如何有效、快速和准确的对点云数据进行分类成了众多学者研究的重点。近年来，针对三维激光点云数据的分类问题，国内外众多学者展开了相关研究。例如：李真等[3]利用K-means算法对不同种类的树进行试验，通过计算高度和强度直方图确定K值和初始聚类中心点来做聚类分析，实现了单株树木原始点云数据分类；胡健等[4]利用地面三维激光扫描仪采集雪地的点云数据，采用相关分析获取点云回波强度与入射角度、距离之间的修正函数，最终实现不同类型雪及地物的分类；Zhang et al.[5]基于自然指数函数阈值将点云划分为不同大小的层次簇，并对其形状特征进行提取和编码，实现了城市环境中的地物分类；Lai和Fox[6]先利用Meanshift聚类方法对环境点云进行分割，然后再利用构成物体点云的宏观特征将各物体分类；王书民等[7]提出了结合点云数据的波形信息，利用模糊C均值的分类方法，将研究区域中树木、地面和建筑进行分类；袁夏等[8]提出利用改进的高斯混合模型对地面不同地物进行分类；倪明等[9]以地面拍摄影像为辅助条件，结合高程、纹理等信息进行地面点云分类。根据上述研究可以看出，学者们主要利用K-means算法、Meanshift算法等传统聚类算法对点云进行分类。但是这些算法都需建立在样本空间为凸的条件下进行分类，当样本空间不为凸时会陷入局部最优，导致分类效果不佳[10]。因此，为了能在任意的样本空间中进行聚类且收敛于最优解，学者们开始采用谱聚类算法进行聚类[11]。谱聚类算法的核心思想来源于谱图理论，由于其具有在任意形状的样本空间上进行聚类并收敛于全局最优解等传统聚类算法难以实现的特点，因此被广泛应用于人脸识别、文本分类和计算机视觉等多个领域[12-14]。

笔者基于Python平台，以三维激光点云数据为数据源，利用ARCGIS软件将点云数据转换为栅格数据，并结合纹理以及形状信息构建谱聚类算法模型，对实验区内的典型地物进行分类，与传统的K-means算法和高斯混合模型进行精度对比分析，探讨谱聚类算法用于典型地物的分类可行性。

1 方法

1.1 谱聚类算法原理

谱聚类算法是聚类分析中的一种新型聚类算法，相较于传统的聚类算法，它建立在谱图理论基础之上[15]。其基本思想：首先利用高斯核函数计算数据之间的欧式距离构建相似度矩阵，将聚类问题转换为图分割问题；之后根据相似度矩阵来获取拉普拉斯矩阵，并利用一种连续的松弛形式将图分割问题转化为拉普拉斯矩阵的谱分解问题，即计算拉普拉斯矩阵前K个最小特征值对应的特征向量；最后利用经典的聚类算法(如K-means算法)对其进行聚类。

(1)相似度矩阵、度矩阵、拉普拉斯矩阵

谱聚类算法的相似度矩阵是通过数据之间的相似度来构建的，其步骤：根据给定的数据集X，X={x1,x2,…,xn}，将每个数据点xi看成无向图G=(V,E)的一个顶点vi，并根据vi和vj连成的边Eij赋权值wij，可以得到无向加权图G(V,E,W)。通过该方法可以把数据集X转化为一个带有权值的无向图G(V,E,W)，并得到相似度矩阵W。相似度是利用高斯核函数来计算的，公式为：

(1)

式中：xi表示第i个点;‖xi-xj‖2表示为xi、xj之间的欧式距离;σ为高斯核参数，控制样本点的邻域宽度，σ越大表示样本点与距离较远的样本点的相似度越大，反之亦然。

将相似度矩阵W的每行元素相加可以得到度矩阵D，公式为：

(2)

式中：D为di组成的n×n对角矩阵，D=diag(d1,d2, …,dn)。

虽然相似度矩阵包含了聚类所需要的数据信息，然而在实际过程中，图分割的过程都是通过求解拉普拉斯矩阵的特征值与特征向量来完成的。拉普拉斯矩阵L可用下式求得：

L=D-W

(3)

(2)图的分割

通过上述分析，利用高斯核函数构造的相似度矩阵完成了从聚类到图分割的过程。那么，怎么对图进行最优分割是谱聚类算法的一个难点。学者们研究出比较好的处理方法是将图分割问题转化为求解连续的拉普拉斯矩阵特征向量问题。即根据拉普拉斯矩阵求解的前K个最小特征值对应的特征向量，可以构成相应的向量空间，其中不同的向量空间对应着不同的分割，即可以达到图分割的目的。目前常见的图分割准则主要有最小割集准则、比例割集准则、规范割集准则和平均割集准则等。本文所使用的图分割准则为比例割集准则，原理：

令图G=(V，E)，图G可分为A、B两部分，并有A∩B=φ,A∪B=V，两个节点间边的权值为W，那么图G划分为A、B两部分的代价函数为：

(4)

其中|A|、|B|表示子图A、B中节点的个数。

(3)算法步骤

①根据给定的数据集X，X={x1,x2,…,xm}，利用高斯核函数计算相似度矩阵W，并构造拉普拉斯矩阵L；

②求取拉普拉斯矩阵L的前K个最小特征值，并计算对应的特征向量u1,u2,…,uk，构成矩阵U={u1,u2,…,uk}；

③取yi为矩阵U的第i行向量，并利用其他经典聚类算法(如K-means算法)，将新样本点Y={y1,y2,…,yn}聚类成簇C1,C2,…,Ck。

1.2 谱聚类算法建模过程

笔者利用谱聚类算法处理高维空间数据的优势，将其应用到三维激光点云数据中，构建基于谱聚类算法的点云数据分类方法。具体建模过程：利用sklearn模块(scikit-learn是python中的第三方库，库中包括分类、回归、聚类和降维等相关算法)中的谱聚类代码集成到遥感影像谱聚类算法中，分别将仅含反射强度特征和加入了纹理、形状信息的多源特征建立谱聚类模型，利用混淆矩阵来评价谱聚类算法模型的分类结果，并根据总体分类精度和Kappa系数来判断其分类效果(图1)。

2 数据采集与处理

2.1 数据采集

以吉林省长春市某一场地为实验区，实验区面积为100 m×100 m，利用Z+F IMAGER 5010C扫描仪对实验区进行数据采集，采集方式为多个测站进行独立扫描，共扫描了14站，实验区包括建筑用地、林地、裸地和草地等地物(图2)。

2.2 特征分析

利用获取的实验区点云数据，导出地物的反射强度值，并随机选取同一地物样本各20个(表1)，分析不同地物反射强度值差异性。从图3中可以看出，不同地物的反射强度值不同且其反射强度值都在一定范围上下波动，验证了反射强度可以作为分类信息将草地、建筑用地、林地和裸地等不同地物区分。

2.3 点云栅格化

在特征分析过程中得出了反射强度可以作为分类信息区分草地、建筑用地、林地和裸地等不同地物，故采用点云数据中的反射强度信息转换为栅格数据。因此利用栅格数据来做后续研究，将点云数据典型地物的分类问题转换为栅格数据的分类问题，利用相关算法对栅格影像进行分类研究。图4为经过拼接、去噪后的点云数据，图5为转换后的栅格数据。

图2 实验区示意图Fig.2 Experimental area diagram

表1 典型地物反射强度值

2.4 特征信息提取

(1)纹理信息提取

在利用栅格影像进行分类时，由于物体的边缘呈现灰度的不连续性和混合像元的存在，仅依靠反射强度信息通常存在丢失像元的现象，而将纹理和反射强度信息结合起来进行分析可以更好地体现地物的宏观性质和细微信息。本文采用统计方法中的灰度共生矩阵(gray level co-occurrence matrix, GLCM) 方式进行分析。GLCM的纹理分析方法由西班牙数学家Haralick于1973年首先提出[16]，GLCM为方阵，维数即灰度级数，通过计算各灰度级之间的联合条件概率密度P(i,j|d,θ)来表示纹理，公式：

图3 实验区地物反射强度折线图Fig.3 Line graph of ground object reflectance in experimental area

图4 点云数据Fig.4 Point cloud data

图5 栅格数据Fig.5 Raster data

P(i,j|d,θ)=φ{(x,y)|f(x,y)=i,f(x+dx,y+dy)=j}

x,y=0, 1, …,N-1

(5)

本文选用的是GLCM纹理分析方法中的方差(variance)和平均值(mean)来进行图像纹理特征提取，提取结果如图6、7所示。

(2)形状信息提取

在栅格影像中，地物的内部信息丰富，使得在统计反射强度和纹理信息时，由于细节过多会出现大量的“盐分”现象，进而导致图像中出现许多空白点和大量的噪声点。本次实验采用高斯滤波方法对影像进行滤波处理并生成0～255的灰度图像，通过直方图分析设置阈值进行边缘信息提取，生成边缘密度二值图像(图8)。

图6 方差Fig.6 Variance

图7 平均值Fig.7 Mean

图8 形状信息Fig.8 Shape information

为了实验更充分，将提取的纹理、形状和反射强度信息进行合成，分别对只包含反射强度信息和多源信息融合后的图像进行分类，对比分类结果。

2.5 分类结果

本文利用谱聚类算法对实验区数据进行分类，为了对比，分别运用了K-means算法和高斯混合模型。试验中为了测试边缘形状、图像纹理等信息对点云数据分类的影响，分别加入了形状信息和纹理信息(图9)。

对比不同分类结果图可以看出，传统的K-means算法和高斯混合模型都很难分出裸地和草地两种相近的地物，并且存在漏分和错分现象，而谱聚类算法表现的更加优越，对于线性地物、不同纹理的分界线都实现了较高的识别；同时，对于反射强度和纹理相近的地物类型，如林地和草地，谱聚类算法分类的效果更好，漏分、错分现象明显减少，谱聚类对数据分布的适应性更强，聚类效果更好。

a.K-means(反射强度信息);b.高斯混合(反射强度信息);c.谱聚类(反射强度信息);d.K-means(多源信息); e.高斯混合(多源信息);(f)谱聚类(多源信息)。图9 分类结果对比图Fig.9 Comparison of classification results

3 精度验证

本次实验进行分类结果精度评价的方法是混淆矩阵[17]，在图像上对应不同地物随机选取1 682个样本进行实验，对只含有反射强度信息和加入了纹理信息、形状信息的分类结果分别计算混淆矩阵中的用户精度、总体分类精度和Kappa系数等指标(表2、图10)。

表2 分类精度指标对比表

图10 分类精度对比图Fig.10 Comparison chart of classification accuracy

从表2可以看出：谱聚类算法对不同地物的分类精度均高于K-means算法和高斯混合模型，其中谱聚类算法对只含有反射强度信息的分类结果的总体分类精度为80.61%，Kappa系数为0.698 2，在结合纹理和形状信息后总体分类精度和Kappa系数分别达到81.36%和0.713 8，验证了谱聚类算法用于三维激光点云数据分类的高效性和适用性，且结合形状和纹理信息后的多源信息图像分类精度会高于仅含反射强度信息数据的分类精度；从地物类型来看，建筑用地和草地的分类精度较高，林地和裸地的分类精度一般，通过结合形状信息和纹理信息后分类精度会有一定程度的提高，其中总体分类精度和Kappa系数分别提高了0.75%和1.56。综上，基于遥感图像的地物分类算法同样适用于点云数据转换成的栅格影像，并且谱聚类算法对点云数据栅格影像的分类精度较高，可行性和适用性都较高，该方法为典型地物分类提供了一个新途径。

4 结论

(1)对于实验区内的草地、建筑用地、林地和裸地等地物，基于反射强度、纹理和形状的多源信息相较于仅基于反射强度信息的点云分类，可以更有效地进行地类识别。

(2)谱聚类算法对点云数据的分类效果较好，其总体分类精度和Kappa系数都高于传统的K-means算法和高斯混合模型。但谱聚类算法耗时较长，后续的研究需对其参数进行一定的优化。

(3)将点云数据转换为栅格数据后再进行分类是一种可行的方法，对地面激光点云数据分类有一定的参考价值。