CN111259944A - 基于快速PCA算法和K-means聚类算法的块石形状分类方法 - Google Patents

Abstract

本发明涉及基于快速PCA算法和K‑means聚类算法的块石形状分类方法，包括以下步骤：(1)根据通过数字图像处理技术获得的只有轮廓信息土石混合体断面CAD图，统计出块石所有边缘轮廓点的坐标及重心坐标；(2)对块石的位置、大小、朝向做统一化处理；(4)采用快速PCA算法对统一化处理所得的块石轮廓半径统计表的每个块石的360条半径为行向量，将其降至二维；(5)选用K‑means聚类算法对降维后的矩阵进行聚类，从而把形状特征相似的块石归为一类；(6)根据定义的形状特征参数对部分分类后的结果进行微调，把具有共性的块石归为一类。该方法能够快速、准确地对大量块石进行移动、旋转、放缩、变向等组合操作，实现大量形状各异的块石的高效分类。

Description

基于快速PCA算法和K-means聚类算法的块石形状分类方法

技术领域

本发明涉及岩土工程数值模拟的技术领域，具体涉及一种基于快速PCA算法和K-means聚类算法的块石形状分类方法。

背景技术

土石混合体是自然界广泛存在的第四纪堆积/沉积物，也是诸多岩土工程的填筑材料，由尺度相对较大、强度较高的块石和相对较细颗粒的弱“土体”组成。由于块石尺度较大，土石混合体的物理力学性质不易采用室内小尺度试验测得，也很难通过开展现场大型试验获得。因此，数值模拟已成为研究土石混合体力学性质的重要手段，且建模的精细程度决定了模拟结果的准确性。故建立考虑块石形状统计规律的数值模型是解决这类问题的基础和关键。

然而，现有的块石分类方法大多只考虑到粒径因素，或者以简单多边形为基准划分为几类，但是由于实际工程是不同地区、不同地质条件下的土石混合体，如果只是简单的将所有的块石轮廓视为简单的多边形及其组合，就可能会导致模拟的结果与实际差异较大。

现有的块石分类方法对于形状复杂、成分不均衡、难以观测的块石应用受限，很难反映土石混合体内部块石之间和块石与细粒土体之间的相互作用，以及块石的扁平度、尖锐度的多因素的共同作用。

发明内容

有鉴于此，本发明提供一种基于快速PCA算法和K-means聚类算法的块石形状分类方法，能够解决现有的土石混合体对于块石分类方面研究不充分的问题。

为此，本发明提出一种基于快速PCA算法的块石形状分类方法，包括以下步骤：

(1)根据通过数字图像处理技术获得的只有轮廓信息土石混合体断面CAD图，统计出块石所有边缘轮廓点的坐标及重心坐标；

(2)采用块石边缘轮廓形状特征统一化程序对块石大小(半径)、形心坐标信息进行处理，将所有块石的重心移到坐标原点(0,0)，确定统一化后的块石的位置及大小；

(3)将步骤(2)获得的位置及大小统一化后的每个块石图形进行旋转处理，使其最长轴落到x轴上，然后对块石轮廓的朝向做统一化处理，并将块石边缘轮廓半径按照同样的顺序输出(这里的同样的意思是指所有块石都按照从0°顺时针选择到-360°的顺序(块石与块石之间的输出顺序同样)，以此来保证相同角度对应每块块石的此角度下的半径，此时将所有块石统一为重心在坐标原点，过重心的长轴与x轴重合，过重心的短轴与y轴重合，过重心的最长半径为1，较尖锐的一方朝向x轴正方向，并将半径按照统一的顺序输出为单行向量，得到以每个块石的360条轮廓半径为行向量的多个块石构成的块石轮廓半径统计表；

(4)采用快速PCA算法对统一化处理所得的块石轮廓半径统计表，即所有块石的半径分布的表格(本实施例中有18块块石，则表格为360*18)的每个块石的360条半径为行向量，将其降维至由只有两个元素构成的二维行向量；

(5)选用K-means聚类算法对降维后的矩阵进行聚类，从而把形状特征相似的块石归为一类，得到聚类结果表。

(6)以某个形状特征参数为基准，把聚类后的同类石块中差异过大的石块挑出来单独归类，把具有共性的块石归为一类，获得最终的块石边缘轮廓半径分类表，即完成块石形状分类。根据定义的形状特征参数，定义在具体实施方式的如图12、如图13所示中有具体说明，指的是扁平度和尖锐度。部分结果是指分类后可能会出现个别块石与其他块石差异巨大并且被单独归为一类的情况，由于本申请的目的是在数量庞大的块石中找出具有共性的块石类群，因此对于这种极其特殊的块石予以舍弃，将他们归为其他类。微调主要是指两个方面，一是如上句所说的，舍弃自己单独属于一类的极其特殊的块石，二是通过计算出的形状特征参数，对与同一类别的其他块石差异过大的块石也进行剔除。对部分分类后的结果进行微调，把具有共性的块石归为一类。

进一步，所述步骤(1)具体为：

块石图形边缘轮廓坐标点提取，同时打开CAD和Excel，加载CAD插件，选取每块块石的边缘轮廓坐标，以每列作为一个完整块石的数据，每完成一块块石坐标的提取，就换到新的列。继续提取，直至完成所有块石坐标点的提取。块石图形是经过现有数字图像技术处理后得到的图形，利用插件统计坐标能明显提高工作效率，通过这两个软件的目的是把大量的块石图形转换成点的坐标的形式，实现从图形到数据的转换，以数据作为操作对象。

进一步，所述步骤(2)具体为：

(21)通过MATLAB编程，对从步骤(1)中提取的坐标进行半径归一化处理，首先将块石所有的x、y坐标减去重心的坐标x_r、y_r，将所有块石的重心移到坐标原点(0,0)，然后将处理后的坐标x’、y’，除以x’的绝对值的最大值与y’的绝对值的最大值中的较大值R_m，保证块石轮廓所有的横纵坐标都小于1。

(22)把步骤(21)处理后的块石视为多条形段组成的多边形，以R＝1.5R_m，重心为圆心，作能够完全包围块石图形的圆，θ∈(0,360°)，每隔1°，画出一条过圆心(相应石块的重心)的射线，让射线与多边形的每条边的线段联立方程组求解交点，得到块石边缘轮廓的360条轮廓半径。

具体包括以下内容：

块石边缘轮廓形状特征统一化程序，利用编写的MATLAB程序，将步骤(1)里面的块石边缘轮廓坐标分别减去对应图形重心的坐标，使重心平移到坐标原点，实现位置(形心坐标)的统一化处理。再把这些坐标除以对应块石的x’坐标的绝对值的最大值与y’坐标的绝对值的最大值中的较大值R_m，使其保证图形的整体形状不变，并实现大小的统一化处理。

用MATLAB编写的以1°为分割单位的半径提取程序由于步骤(21)的图形半径已经经过归一化处理，因此能够完全包围块石图形的圆的半径为

这里取1.5R_m，通过直线方程求交点坐标的方法来求出圆和块石图形的交点坐标，确定统一化后的块石半径大小。特例，当块石图形的边界线与Y轴平行的时候，将相邻点的X坐标分别±0.0001，保证直线方程有解。

进一步，所述步骤(3)具体为：

(31)将步骤(22)所得的块石边缘轮廓的360条半径。按照与坐标原点三点共线的原则分成180组，求出每组两个交点的距离。

(32)找出距离最大点所在的位置M，并将图形以重心为旋转中心，过重心的最长轴与x轴正方向所成的角为旋转角m，顺时针旋转m°，其角度与位置相对应的作用。

(33)求出旋转后的每个点与坐标原点所连接成的直线的反正切值，并将其划分为正负两组，以1°的圆心角分割为360个点围成的区域，分别按照0～-180°和180°～0的顺序进行排序，使每个点与坐标原点所连接成的直线的反正切值按照0°顺时针旋转到-359°的顺序排序。

这是由于反正切值函数y＝arctan x的值域是(-π/2，π/2)，并且单调区间为(-∞，+∞)递增，故而，将算出来的每个点与坐标原点所连接成的直线的反正切值分成两个部分，一部分的y值从(0，-π/2)按递减的顺序排序，其对应的角度也必然是顺时针从0旋转到-180°；另一部分的y值从(-π/2，-π)按递减的顺序排序，其对应的角度也必然是顺时针从-180°旋转到-359°，这一步的意义因为在后续的块石边缘轮廓图形生成过程中，本申请所采取的顺序同样是从顺时针从0旋转到-359°，这里体现了前后的相对应的特点。

(34)找出经过(33)处理后的图形的块石轮廓的上、下、左、右四个参考点。分别求出上下两点与左、右两个端点的夹角并进行比较，如果右夹角的度数大于左夹角，则让图形关于y轴对称，以此来保证较尖锐的角始终朝向x轴正方向，也就是变向。如果需要变向(变向是指保证右边的夹角度数小于左边的夹角度数，如果右边的夹角度数大于左边的夹角度数，也就是右边没有左边尖锐，那么才进行变向；反之，则不变)，返回参考(33)用同样的方法和顺序进行排序，以保证结果的准确性；如果不需要变向则保持(33)的结果不变，把旋转朝向统一后的块石边缘轮廓半径统计表保存到Excel表格中。

具体包括以下内容：

块石最长轴查找程序(步骤31)，将步骤(22)所得的块石边缘轮廓的360份细化坐标按照与坐标原点三点共线的原则分成180组，得到了180条过重心的直径，利用两点距离公式求出所有的直径。

块石轮廓旋转程序(步骤32)，用MATLAB内置的find函数找到最大值所在的位置M，通过坐标旋转公式，将所有的坐标顺时针旋转m°，最长轴顺时针旋转m°就能使最长轴落在X轴上，从而实现块石角度的统一化处理。

块石方向改变程序(步骤34)，找出旋转后的块石轮廓的对应角度为0°、180°的点坐标，以及单个块石轮廓矩阵y坐标最大及y坐标最小的点坐标，这样就分别得到了右、左、上、下四个参考点。分别求出右点与上下两点连线所成的角(右夹角)以及左点与上下两点连线所成的角(左夹角)，比较两个角的大小，如果右夹角的度数大于左夹角，则将对应图形的所有横坐标乘以-1，使其关于y轴对称。

进一步，所述步骤(4)具体为：

(41)利用快速PCA算法对块石轮廓半径进行降维，从而将块石轮廓半径数据从高维空间投影到低维空间中，并尽可能的在低维空间中表示原始半径关系。

(42)画出降低到二维后的块石边缘轮廓半径关系图(即降维后的各个行向量的关系)，并初步判断和预估划分的种类数(K值)。

具体包括以下内容：

快速PCA算法的主要步骤为，以每个块石轮廓半径分布为行向量，计算每个行向量的均值和轮廓矩阵的协方差矩阵，然后计算协方差矩阵转置的前两个特征值和特征向量，将多维降成二维，然后用归一化处理后得到块石轮廓矩阵乘以特征向量，并把特征向量归一化为单位特征向量就可以实现降维。

进一步，所述步骤(5)具体为：

通过MATLAB的plot函数，以二维平面直角坐标图的方式表示出降维后所有坐标点的关系，进行初步观察判断，估计需要聚类的数目K，原则是使得同类的元素之间的距离最小，而类与类之间的距离最大。

(51)通过MATLAB程序画出层次聚类图，并通过进一步的观察缩小K值的范围(通过观察缩小K值，看画出来的层次聚类图上面点的距离，距离越近的几个点在一类的可能性就越大)，其中每从垂直于y轴的方向做一条直线与聚类树的交点为一类。

(52)继续通过MATLAB程序作出K-SSE(聚类—每类内部所有块石到聚类中心的距离平方和)曲线，利用肘部法则找到K值畸变程度的改善效果下降幅度最大的位置，确定最终的K值。

(53)输入步骤(52)获得的K值，通过K-means聚类算法进行聚类，并将输出结果以块石边缘轮廓的分类图及聚类结果表的形式表示出来便于查看。

前K个的选取方法主要是根据(42)中所述的降低到二维后的块石轮廓半径关系图以及来确定，这里是一个初步确定，在后续步骤(5)中，多次选取的目的是缩小K的范围，提高了分类数量的准确度，并且用多种方式选取的目的是对单一模式加以验证。此外，对K值如此小心翼翼多次确定的原因是K-means算法的准确度非常依赖K值。

具体包括以下内容：

层次聚类程序(层次聚类程序为现有的程序，用来让K值的确定更加准确，整个过程并不会直接利用层次聚类的任何结果值)，首先，把块石轮廓矩阵的每行都单独看做是一类，其次，计算出所有行向量两两之间的距离，构成距离矩阵，并把距离最近的两类合并为一个新类，然后，计算出新类与当前各类的距离，再不断的进行行向量的合并与距离计算，直到只有一类为止。

K-means算法的参数是每一类的中心位置和其内部观测值的位置，其目标是是以成本函数最小化。K-means成本函数公式如下：

其中μ₁,μ₂,…,μ_ik(ik＝1,2,…,K)来表示聚类中心，用c(1),c(2),...,c(jk)来存储与第jk个块石轮廓矩阵行向量最近的聚类中心的索引，μ_ik是第ik个类的中心位置，即平均值，成本函数是各个类畸变程度之和，ik表示的是类的数量，jk表示的是类里面块石的数量，ik是一个大类，jk是ik里面的分支，式中的竖线数学意义是表示两个向量的距离，即图7的两个点的距离。

可以通过作出K-SSE(聚类—每类内部所有块石到聚类中心的距离平方和)曲线，用肘部法则来估计聚类数量，随着k值的增大，平均畸变程度会减小；每个类包含的块石数量会减少，于是每个块石距离中心会更近，畸变程度的改善效果下降幅度最大的位置所对应的k值称为肘部。

上述K值的确定可以作为参考，提倡的做法是从实际问题出发，人工指定比较合理的K值，通过多次随机初始化聚类中心选取比较满意的结果。

K-means是一个迭代算法，在降维后的块石轮廓矩阵中，首先选择K个随机的点，称为聚类中心；以块石轮廓矩阵中的每行作为行向量，按照距离K个中心点的距离，将其与距离最近的中心点关联起来，与同一个中心点关联的所有块石聚成一类。计算每一个类的平均值，并将该类所关联的中心点移动到平均值的位置。重复上述步骤，直至中心点不再变化。

进一步，所述步骤(6)具体为：

(61)根据保存的聚类结果表和块石边缘轮廓半径统计表，以及块石边缘轮廓分类图，求出所有图形的几个形状特征参数，获得相应的形状特征参数表格，如扁平度和尖锐度等。

(62)结合得到的形状特征参数表格和降维后的块石边缘轮廓半径分布图，把与同簇其他块石差异过大的块石挑选出来单独归类，得到微调后的块石边缘轮廓半径分类表。

(63)根据微调(微调指的是经(62)步骤处理后得到的)后的块石边缘轮廓半径分类表，把具有共性的块石归为一类，把与其他几簇差异都比较大的单独罗列，画出块石的分类结果图便于建模时选用。

具体包括以下内容：

块石形状特征参数计算程序，主要计算块石的扁平度和块石尖锐度，计算公式如下：

其中l_max,b_max分别为块石经步骤(3)统一化处理后的过重心的最长轴、最短轴(块石图形最上与最下两个点的纵坐标之差)。

其中θ_a、θ_b、θ_c、θ_d分别为块石图形的右、左、下、上四个点连线形成的夹角的角度值；α₁、α₂分别代表水平方向尖锐度，竖直方向尖锐度。

块石图形生成程序，根据块石的形状特征参数表、块石的形状分类图、降维后的块石边缘轮廓半径分布图这几个图表，对分组结果进行微调，并画出新的分类图。

与现有技术相比，本发明的有益效果为：

本发明基于计算机机器学习领域的快速PCA降维算法和K-means聚类算法，对块石的边缘轮廓形态特征进行快速、准确的分类，减少了人为干预、运行效率高。

本发明方法定义了块石的部分形状特征参数并进行微调，避免了部分与簇内和其他簇的岩石都不具有共性而被错误归为一类的误差。

本发明为了对块石形状进行分类，需要对块石统一处理，减少干扰因素的影响，通过对每个块石进行移心、旋转、归一、变向等标准化处理，使其统一为重心在坐标原点，过重心的长、短轴(最长、最短直径)与x轴、y轴重合，过重心的最长半径为1，较尖锐的一方朝向x轴正方向，并将半径按照统一的顺序输出为单行向量。本发明加入对于尖锐角朝向的处理，其目的是避免两个块石轮廓相似却在旋转过程中发生的朝向改变被错认为不同形状的影响。基于MATLAB，编写了块石形状分类处理软件，只需要导入块石的块石边缘轮廓特征点坐标，便可快速进行形状分类，并且可以导出处理后的结果作用于后续的建模，操作简单、便捷、拓展性强。

通过上述处理，可将块石的二维轮廓转化为半径ρ，角度θ的函数，鉴于角度都是按照0°～359°的顺时针顺序，可以只比较半径ρ方向的相似程度，考虑到人工方法费时、费力且具有较强的主观性，采用K-means聚类算法进行聚类分类。

K-means聚类算法具有计算速度快、可拓展性好的优点，其步骤是，随机选取K个初始点作为聚类中心(初次选取的聚类中心称为种子聚类中心)，然后计算每个块石与各个种子聚类中心之间的距离，把每个块石分配给距离它最近的聚类中心，聚类中心以及分配给它们的块石就代表一个聚类。每分配一个块石样本，聚类的聚类中心会根据聚类中现有的块石样本重新计算。这个过程将不断重复直到满足已经设定的终止条件。终止条件可以是没有(或最小数目)的块石对象被重新分配给不同的聚类，或者没有(或最小数目)聚类中心再发生变化。

本申请先将轮廓矩阵的维度降低，然后再进行聚类，使数据可视化且降低计算机的计算负担，快速PCA算法可以从多维块石轮廓矩阵中解析出主要影响因素，对于维度较高的块石轮廓矩阵，用块石轮廓矩阵的转置矩阵的协方差矩阵来代替协方差矩阵，将快速PCA降维算法与K-means聚类算法结合起来应用于土石混合体块石的分类，减小了计算量，提高了计算速度。

本申请根据聚类分类的结果，对块石的形状特征参数加以总结，找出块石的扁平度尖锐度都差距较大的结果，将具有相似特征的块石图形归为一类，与每一类差别都较大的块石图形挑拣出去，作为其他类。

本发明方法对通过数字图像处理技术得到的真实块石形状进行研究，保证了后续块石形态分类更具有真实性、可靠性。

附图说明

图1为本发明的方法流程图

图2为本发明的通过数值图像处理获得的待分类的18个块石图

图3为本发明的以坐标点表示的某个块石图形分割示意图

图4为本发明的块石轮廓半径与角度关系图

图5为本发明的块石的重心与最长轴图

图6为本发明的旋转前后的块石轮廓对比图

图7为本发明的快速PCA算法降至二维的向量关系图

图8为本发明的土石混合体集合层次聚类图

图9为本发明的K-SSE曲线图

图10为本发明的块石聚类时不同K值的轮廓系数图

图11为本发明的程序聚类结果图

图12为本发明的块石图形边缘轮廓尖锐度示意图

图13为本发明的微调后的聚类结果图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明基于快速PCA算法和K-means聚类算法的块石形状分类方法流程图(参见图1)，分为三个部分，包括：

(1)获取块石形状信息

首先，对通过数字图像技术处理得到的二维土石混合体块石边缘轮廓MATLAB图，块石图形随机分布图如图2所示，并将其导入到CAD中使用CAD插件，提取出块石图形每个块石边缘轮廓点的坐标；然后根据块石图形边缘轮廓半径提取程序，对块石图形边缘轮廓点进行分割，使其按照1°的单位分割为360份(如图3为单个图形的分割示例)，并将其按照同样的顺序输出，最后保存处理后的360等分块石边缘轮廓坐标以及半径。

(2)块石形状统一处理

首先，要实现位置(形心)的统一化处理，将块石图形特征坐标分别减去对应图形重心的坐标，使重心平移到坐标原点；

其次，要实现大小的统一化处理，将重心平移后的块石图形坐标除以对应重心平移后块石的横、纵坐标的绝对值的最大值的较大值R_m，使其保证图形的整体形状不变；

然后，要实现倾角的统一化处理，将所有过重心的半径两两对应分成180组，利用两点距离公式求出最大的直径，再用MATLAB内置的find函数找到最大值所在的位置M，过重心的最长轴与x轴正方向所成的角为旋转角m，则将最长轴顺时针旋转m°就能使最长轴落在X轴上(如图5、图6所示)；

最后，要实现朝向的统一化处理，找出旋转后的块石轮廓的上、下、左、右四个参考点。分别求出上下两点与左、右两个端点的夹角并进行比较，如果右夹角的度数大于左夹角，则让图形关于y轴对称，以此来保证较尖锐的角始终朝向x轴正方向，此时所有块石统一为重心在坐标原点，过重心的长轴与x轴重合，过重心的短轴与y轴重合，过重心的最长半径为1，较尖锐的一方朝向x轴正方向，并将半径按照统一的顺序输出为单行向量，得到以每个块石的360条轮廓半径为行向量的统一处理后的多个块石构成的块石边缘轮廓半径表，记为表2。

(3)块石分类结果输出

首先，将块石边缘轮廓半径表(表2)以每行为一个类别，其中每行都包含了一个完整的块石半径向量，通过快速PCA主成分分析，对每个行向量都实施降维，使其降低到二维，得到降维后的块石边缘轮廓半径表，记为表6，便于后续的计算处理，并且还能够实现数据的可视化(参见图7)，得到降维后的块石边缘轮廓半径分布图，图7中X、Y分别代表经过PCA主成分分析的第一主成分、第二主成分，它代表了由原始数据中的变量经过线性变换得到的新变量；图中圆圈的数字，表示的是与X，Y相对应的块石的编号，两个点距离越近表示对应块石的相似性越强。

然后，通过编写的MATLAB程序画出层次聚类图(图8中以欧式距离为纵坐标，以数据集中的元素为纵坐标画出数据集的18(块石总数，也就是块石顺序编号)个元素层次聚类竖图。数据集是指降维后的块石边缘轮廓半径表(表6)的数据。

和K-SSE(聚类—每类内部所有块石到聚类中心的距离平方和)曲线(参见图9)，进一步缩小K的取值范围，也可以借助图10所示的不同类别对应的块石平均轮廓系数图，进一步缩小K的取值范围；上述K值的确定可以作为参考，最为提倡的做法是根据实际问题来确定比较合理的K值，并通过多次随机初始化聚类中心选取比较满意的结果；最后，通过K-means聚类算法将边缘轮廓比较相似的矩阵聚成一类，并将输出结果以图形(图11)的形式表示出来便于查看，同时程序里也包含了以表格的形式表达分类结果的命令，即获得聚类结果表，便于进一步完善。

如图2所示，用数字图像处理得到的土石混合体部分块石图形随机分布图，将处理后的块石图形导入到CAD中，并且打开CAD和Excel，加载CAD插件并利用CAD插件对块石图形边缘轮廓的交点进行提取，然后将土石混合体的半径和坐标信息导出到Excel，生成原始块石边缘轮廓半径表，记为表1(本申请中未给出)，便于后续的使用和处理。

包括以下内容：

1)利用数字图像处理技术将块石图片二值化并提取块石边界线。

2)将处理后的块石图形导入到CAD中，并利用CAD插件提取块石边缘轮廓的交点坐标。

将处理后得到的块石边缘轮廓图形坐标保存成CSV格式的表格，记为表1，便于后续的使用。

如图3所示，任取一块块石的轮廓提取图，其中对每个块石图形的边界分别进行分割处理，以1°为单位，每个块石提取出360条半径。

包括以下内容：

求出每个块石图形的重心，由于块石图形的形状不规则，既有凸多边形又有凹多边形，因此通过公式(1)、公式(2)来求出重心。

其中，C_ix、C_iy表示图形被分割的每个小三角形的重心横、纵坐标；A_i表示这些小三角形的面积，n＝360。

求图形的半径分布和点的坐标，首先将所有的x、y坐标减去重心的坐标x_r、y_r；然后将处理后的坐标x’、y’，除以所有x’、y’的绝对值的最大值的较大值R_m，经过处理后能够完全包围块石图形的圆的半径为

(取1.5R_m)。以R＝1.5R_m，，θ∈(0,360°)，每隔1°，画出一条过圆心的射线与块石边缘轮廓相交，通过求两条直线方程交点坐标的方法来求出块石边缘轮廓分割点。

图4是一个完整的块石图形的块石轮廓半径与角度的关系图，x轴是角度按照顺时针1°到-360°，y轴是与角度对应的半径，由于步骤2中已经得到轮廓折线交点的坐标，但是由于块石的形状普遍比较复杂，且不同图形的边缘轮廓交点坐标的数量不同。因此，图4为统一化处理后的半径与角度的平面直角坐标系关系图，由以1°为单位的360个坐标点组成。

包括以下内容：

块石几何特征统一处理程序，利用编写的MATLAB程序，将块石边缘轮廓坐标分别减去对应图形重心的坐标，使重心平移到坐标原点，以将块石进行统一的位置移动，使其重心坐落于原点，实现位置的统一化处理。再把这些坐标除以对应块石的x、y坐标的绝对值的最大值的较大值R_m，使其保证图形的整体形状不变，并实现大小的统一化处理。

用MATLAB编写的以1°为分割单位的半径提取程序，由于块石半径已经经过归一化处理，因此能够完全包围块石图形的圆的半径为

这里取1.5R_m，通过直线方程求交点坐标的方法来求出圆和块石图形的交点坐标，用以表示不同块石直接的距离关系。特例，当块石图形的边界线与Y轴平行的时候，将相邻点的X坐标分别±0.0001，保证直线方程有解。

如图5所示，块石的重心与过重心的最长轴图，求出块石的几何中心(重心)然后找出过重心的最长轴线。

包括以下内容：

将得到的块石边缘轮廓的360份细化坐标按照与坐标原点三点共线的原则分成180组，使用两点距离公式，求出每组两个交点的距离。

找出距离最大点所在的位置M，画出过M和M+180°的直径，找出的最长轴部分。

如图6所示，单个块石图形旋转前后的对比图，将块石图形以重心为旋转中心，过重心的最长直径与x轴正方向所成的夹角为旋转角旋转。

包括以下内容：

用MATLAB内置的find函数找到最大值所在的位置M，求出旋转后的每个点与坐标原点所连接成的直线的反正切值，并将其划分为正负两组，分别按照0～-180°和180°～0的顺序进行排序，并组合起来，则将最长轴逆时针旋转M°就能使最长轴落在X轴上，从而实现块石角度的统一化处理。

找出旋转后的块石轮廓的对应角度为0°、180°的点坐标，以及单个块石轮廓矩阵y坐标最大及y坐标最小的点坐标，这样就分别得到了右、左、上、下四个参考点。分别求出右点与上下两点连线所成的角(右夹角)以及左点与与上下两点连线所成的角(左夹角)，比较两个角的大小，如果右夹角的度数大于左夹角，则将对应图形的所有横坐标乘以-1，使其关于y轴对称，从而实现朝向的统一化处理。

如图7所示，基于PCA降维技术，根据输入的统一化处理后的18个块石边缘轮廓半径统计表，记为表2，重新找到更能描述各个点之间的距离关系的正交的坐标轴，通过fast-PCA函数来对块石边缘轮廓矩阵进行快速主成分分析和降低到二维处理，输出块石边缘轮廓矩阵特征向量组成的降维后的矩阵，其中每个行向量一个块石边缘轮廓样本，列数为降维后的块石边缘轮廓特征维数(也就是2)。

包括以下内容：

快速PCA算法的主要步骤为，以每个块石轮廓半径分布为行向量，计算每个行向量的均值和轮廓矩阵的协方差矩阵，然后计算协方差矩阵转置的前2个特征值和特征向量，然后用原数据乘以特征向量，并把特征向量归一化为单位特征向量就可以实现降维。

计算公式如下：

假设X_ip(ip＝1,2,…,18)是360个块石边缘轮廓半径组成的行向量的转置

X＝[X₁,X₂,…,X₁₈]^T，并且μ_ip是第ip个(ip＝1,2,…,18)行向量的均值，即μ_ip＝E(X_ip)

则协方差矩阵的第ip,jp项是：

∑_ip,jp＝cov(X_ip,X_jp)＝E[(X_ip-μ_ip)(X_jp-μ_jp)] (3)

通过MATLAB的plot函数，以二维平面直角坐标图的方式表示出降维后所有坐标点的关系(参见图7)，进行初步观察判断，估计需要聚类的数目K，原则是使得同类的元素之间的距离最小，而类与类之间的距离最大。

如图8所示，运用层次聚类算法，计算不同类别数据点间的相似度，创建一棵有层次的嵌套聚类树。其中，不同类别的原始数据点是树的最低层，树的顶层是一个聚类的根节点，相似度的计算采用欧氏距离，见公式(4)、(5)。将同类的数据点进行组合，组合后重新计算距离，例如，当计算集合(C、D)(通过不断的加入和剔除找出类与类之间距离最小，本类与外类距离最大的分配方式)到B的距离时，见公式(6)。

其中x_D、x_B为D、B两点的横坐标，y_D、y_B为D、B两点的纵坐标；DC_B为二者的平均值，其代表C、D集合与B点的距离；D_C为C、B两点之间的距离，D_B为D、B两点之间的距离。

包括以下内容：

层次聚类程序，首先，把经过降维的块石轮廓矩阵的每行都单独看做是一类，其次，计算出所有行向量两两之间的距离，构成距离矩阵，并把距离最近的两类合并为一个新类，然后，计算出新类与当前各类的距离，再不断的进行行向量的合并与距离计算，直到只有一类为止。

如图9所示，K-SSE(聚类—每类内部所有块石到聚类中心的距离平方和)曲线图，用肘部法则(找到斜率突变的地方，术语以人的手肘来形容，所以叫做肘部法则)来估计聚类数量，随着K值的增大，平均畸变程度(y轴的值)会减小；每个类包含的块石数量会减少，于是每个块石距离中心会更近，畸变程度的改善效果下降幅度最大的位置(也就是斜率发生突变的地方)所对应的K值称为肘部。

包括以下内容：

K-means算法的参数是每一类的中心位置和其内部观测值(也就是公式(7)的x_jk)的位置，其目标是是以成本函数(J)最小化。K-means成本函数公式如下：

其中μ₁,μ₂,…,μ_ik(ik＝1,2,…,K)来表示聚类中心，用c(1),c(2),...,c(jk)来存储与第jk个块石轮廓矩阵行向量最近的聚类中心的索引，μ_ik是第ik个类的中心位置，即平均值，成本函数是各个类畸变程度之和，ik表示的是类的数量，jk表示的是类里面块石的数量，ik是一个大类，jk是ik里面的分支，式中的竖线数学意义是表示两个向量的距离，即图7的两个点的距离。

如图10所示，块石的分类数与平均轮廓系数图，轮廓系数结合了聚类的凝聚度和分离度，用于评估聚类的效果。该值处于-1～1之间，值越大，表示聚类效果越好。

包括以下内容：

对于每个块石样本ic，计算降维后的块石半径行向量ic与其同一个簇内(簇a)的所有其他块石半径行向量的距离平均值a(ic)，用于量化簇内的凝聚度。

选取块石样本ic外的另一个簇b，计算ic与b中所有块石半径行向量的平均距离，遍历所有簇，找到最近的这个平均距离,记作b(ic)，用于量化块石不同簇间分离度，计算所有ic的轮廓系数，求出平均值即为当前聚类的整体轮廓系数。

上述几种方法可以作为K值确定的参考，根据实际问题划分一个大致的范围，人工指定比较合理的K值，通过多次随机初始化聚类中心选取比较满意的结果。通过聚类获得聚类结果表即块石边缘轮廓分类表，记为表3。

如图11所示，基于计算机机器学习领域的PCA降维算法和K-means聚类算法的输出结果图，以实现对块石的形态特征进行快速、准确的分类。

包括以下内容：

对聚类后得到的结果矩阵(MATLAB编程是以矩阵为单位，所以初始结果中含有0)进行处理，去除掉里面的零元素。

根据每个矩阵的值与对应的行数列数，用MATLAB在对应的列画出对应的图形(图11)。

为了进一步的归纳和整理，输出块石的分类结果表格(表3)，按照其形状特征参数进行微调。

如图12所示，单块块石的边缘轮廓四个尖锐角的示意图，用以进行块石图形的进一步分类和筛选，找出与簇内其他块石差异大的块石。

包括以下内容：

计算出块石的扁平度(μ_b)，并将计算得到的扁平度块石形状特征参数表保存到Excel中，记为表4a。

计算公式如下：

其中l_max,b_max分别为块石经统一化处理后的过重心的最长轴、短轴(块石图形最上与最下两个点的纵坐标之差)。

计算出块石尖锐度(α₁、α₂)，并将计算得到的尖锐度块石形状特征参数表保存到Excel中，记为表4b中。

计算公式如下：

其中θ_a、θ_b、θ_c、θ_d分别为块石图形的右、左、下、上四个点连线形成的夹角的角度值。扁平度和尖锐度(两者中有一者不满足则需要被剔除)与本类的其他图形差距为20％时应予以剔除，这个数值的选取是根据多次试验选取的较好数值，但是由于可应用块石的种类和范围不同，也可以人为指定这个数值(例如在样本数量少的块石边缘轮廓图形集合中，可以将20％的范围加大，而在样本数量较多的块石边缘轮廓图形集合中，可以将20％的范围缩小)。

如图13所示为微调后的二维块石图形边缘轮廓分类图，综合得到的尖锐度块石形状特征参数表、扁平度块石形状特征参数表和降维后的块石边缘轮廓半径分布图以及块石边缘轮廓分类表(表3)，对分类表进行调整，完成块石形状分类，分类结果见图13。

虽然结合附图描述了本发明的实施方式，但是本领域技术人员可以在不脱离本发明的精神和范围的情况下做出各种修改和变型，这样的修改和变型均落入由所附权利要求所限定的范围之内。

本发明未述及之处适用于现有技术。

Claims (7)

1.一种基于快速PCA算法和K-means聚类算法的块石形状分类方法，其特征在于，包括以下步骤：

(1)根据通过数字图像处理技术获得只有轮廓信息的土石混合体断面CAD图，统计出块石所有边缘轮廓点的坐标及重心坐标；

(2)位置及大小统一化：对块石的形心坐标及半径进行统一化处理，将所有块石的重心移到坐标原点(0,0)，确定统一化后的块石的位置及大小；

(3)旋转朝向统一化：将步骤(2)获得的位置及大小统一化后的每个块石图形进行旋转处理，使其最长轴落到x轴上，然后对块石轮廓的朝向做统一化处理，并将所有块石边缘轮廓半径都按照从0°顺时针选择到-360°的顺序的顺序输出；此时将所有块石统一为重心在坐标原点，过重心的长轴与x轴重合，过重心的短轴与y轴重合，过重心的最长半径为1，较尖锐的一方朝向x轴正方向，并将半径按照统一的顺序输出为单行向量，得到以每个块石的360条轮廓半径为行向量的多个块石构成的块石轮廓半径统计表；

(4)降维：采用快速PCA算法对步骤(3)获得的块石轮廓半径统计表进行降维处理，将其降至二维，得到降维后的块石边缘轮廓半径分布图；

(5)聚类：选用K-means聚类算法对降维后的矩阵进行聚类，从而把形状特征相似的块石归为一类，获得聚类结果表；

(6)确定所有石块的形状特征参数，以某个形状特征参数为基准，把聚类后的同类石块中差异过大的石块挑出来单独归类，把具有共性的块石归为一类，获得最终的块石边缘轮廓半径分类表，即完成块石形状分类。

2.根据权利要求1所述的基于快速PCA算法和K-means聚类算法的块石形状分类方法，其特征在于，在所述步骤(1)中，包括以下步骤：

块石图形边缘轮廓坐标点提取，同时打开CAD和Excel，加载CAD插件，选取每块块石的边缘轮廓坐标，以每列作为一个完整块石的数据，每完成一块块石坐标的提取，就换到新的列继续提取下一块块石坐标，直至完成所有块石坐标点的提取。

3.根据权利要求1所述的基于快速PCA算法和K-means聚类算法的块石形状分类方法，其特征在于，在所述步骤(2)中，包括以下步骤：

(21)利用MATLAB对从步骤(1)中提取的块石坐标进行半径归一化处理，首先将块石所有的x、y坐标均减去重心的坐标x_r、y_r，得到块石处理后的坐标x’、y’，找出所有x’的绝对值的最大值与y’的绝对值的最大值中的较大值，最后将每个处理后的坐标x’、y’，除以该较大值R_m，保证块石轮廓所有的横纵坐标都小于1；

(22)把步骤(21)处理后的块石视为多条形段组成的多边形，以R＝1.5R_m，重心为圆心，作能够完全包围块石图形的圆，θ∈(0,360°)，每隔1°，画出一条过圆心的射线，让射线与多边形的每条边的线段联立方程组求解交点，确定块石边缘轮廓的360条轮廓半径。

4.根据权利要求3所述的基于快速PCA算法和K-means聚类算法的块石形状分类方法，其特征在于，在所述步骤(3)中，包括以下步骤：

(31)确定块石最长轴：将步骤(22)获得的块石边缘轮廓的360条半径按照与重心三点共线的原则分成180组，求出每组两个交点的距离，找出距离最大的两个交点所在位置M，确定为块石的最长轴；

(32)块石轮廓旋转：将块石以重心为旋转中心，过重心的最长轴与x轴正方向所成的角为旋转角m，顺时针旋转m°；

(33)块石排序：求出旋转后的每个点与(0,0)点所连接成的直线的反正切值，并将其划分为正负两组，以1°的圆心角分割为360个点围成的区域，分别按照0～-180°和180°～0的顺序进行排序，使每个点与坐标原点所连接成的直线的反正切值按照0°顺时针旋转到-359°的顺序排序；

(34)块石变向：确定步骤(33)排序后的块石轮廓的上、下、左、右四个参考点，分别求出右点与上下两点连线所成的右夹角以及左点与上下两点连线所成的左夹角，比较左右夹角的大小，若右夹角的度数大于左夹角，则需要变向，使块石图形关于y轴对称，然后返回步骤(33)重新排序，以此来保证较尖锐的角始终朝向x轴正方向；

若右夹角的度数不大于左夹角则不需要变向，保持步骤(33)的结果不变；获得旋转朝向统一后的块石边缘轮廓半径统计表。

5.根据权利要求1所述的基于快速PCA算法和K-means聚类算法的块石形状分类方法，其特征在于，在所述步骤(4)中，包括以下步骤：

(41)利用快速PCA算法对块石轮廓半径进行降维，从而将块石轮廓半径数据从高维空间投影到低维空间中，并尽可能的在低维空间中表示原始半径关系；

(42)画出降低到二维后的块石轮廓半径关系图，并初步判断和预估划分的种类数，即K值。

6.根据权利要求5所述的基于快速PCA算法和K-means聚类算法的块石形状分类方法，其特征在于，在所述步骤(5)中，包括以下步骤：

(51)通过MATLAB画出层次聚类图，并通过观察缩小步骤(42)的K值的范围，其中每从垂直于y轴的方向做一条直线与聚类树的交点为一类；

(52)继续通过MATLAB作出K-SSE曲线，即聚类-每类内部所有块石到聚类中心的距离平方和曲线，利用肘部法则找到K值畸变程度的改善效果下降幅度最大的位置，确定最终的K值；

(53)输入步骤(52)获得的K值，通过K-means聚类算法进行聚类，获得聚类结果表。

7.根据权利要求1所述的基于快速PCA算法和K-means聚类算法的块石形状分类方法，其特征在于，在步骤(6)中所述形状特征参数为扁平度和尖锐度，扁平度μ_b和块石尖锐度，计算公式是：

其中l_max为块石经步骤(3)统一化后的过重心的最长轴，b_max为石经步骤(3)统一化后块石图形最上与最下两个点的纵坐标之差；

其中θ_a、θ_b、θ_c、θ_d分别为块石图形的右、左、下、上四个点连线形成的夹角的角度值，其中θ_a、θ_b为水平方向的两个角，θ_c、θ_d为竖直方向的两个角；α₁、α₂分别代表水平方向尖锐度、竖直方向尖锐度。

Images