一种基于改进K-means的微震信号聚类方法及装置与流程

本发明涉及数据挖掘技术领域，尤其涉及一种基于改进k-means的微震信号聚类方法及装置。

背景技术：

煤矿生产过程中通常存在顶底板断裂、卸压与掘进爆破，均使得岩体受到外力的影响并释放能量而产生持续一段时间的震动。通常对于振动频率小于100hz，事件能量在103～1011j的震动波称为微震。在煤岩变形破裂的过程中所产生的微震信号具有时间短、突变快的特性。矿井爆破和机械振动等外界环境因素的干扰噪声也会在一定程度上影响到微震信号的波形形态，导致信号失真。综上，微震信号的传播过程中会受到很多因素的干扰，故需要对其进行深入研究。

随着开采深度的增加和采空区面积的扩大，矿山微震的发生越来越频繁。微震监测系统在矿山领域中得到了广泛应用，主要包括室内柜式微震监测系统、野外悬挂式微震实时监测系统以及地表浅埋式wifi微震实时监测定位系统。微震监测系统连续不断地产生大量的微震信号，常常也伴随着各种各样的噪声信号。由于微震事件受干扰较多，呈现不同类型，无论是对于微震事件的处理还是对事故发生后复工复产的方案制定都需要对微震事件的类型进行判识。目前国内外研究学者对微震信号的研究已经逐渐深入，通过微震信号的非线性非平稳的特征，小波分析和经验模态分解方法在对微震信号的特征提取研究方面有了很大的进展，微震波到时可以通过滑动平均值法、sta\lta方法、分形方法等方法实现。所提出的这种聚类装置能够实现对微震信号的有效聚类。微震信号聚类分析借助聚类算法的理论支撑，不仅可以发现微震事件的潜在特征还能够有效地降低人工划分的主观性，能够发现微震事件的活动规律和分布特征。由于微震信号是一种非平稳非线性信号，发现其潜在分布特征相对较难，传统的k-means算法仅能够针对数据点集进行分析聚类，因此所提出的装置中将处理时间序列数据的dtw(dynamictimewarping，动态时间归整算法)距离与k-means相结合来对微震事件聚类。综上所述，从大量的微震信号中挖掘有用的信息，对微震事件做出预测性的判断已经成为研究的热点与难点。

技术实现要素：

本发明主要是解决现有技术中所存在的技术问题，从而提供一种适用于微震信号聚类分析、处理效率高、可靠性好的基于改进k-means的微震信号聚类方法及装置。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

本发明提供的基于改进k-means的微震信号聚类方法，其包括以下步骤：

s1、将微震信号传入到k-means控制机中，生成样本数据和初始聚类中心；

s2、计算样本数据与初始聚类中心之间的dtw距离；

s3、比较dtw距离大小，标注标签后进行第一次聚类，将第一次聚类得到的结果传给dba更新器中，得到新的聚类中心；

s4、比较初始聚类中心与更新后的聚类中心是否相等，若是，则执行步骤s5，若否，则执行步骤s6；

s5、输出聚类结果；

s6、返回步骤s3-s4，进行下一次聚类。

进一步地，所述步骤s1之前包括步骤：

s0、对微震信号进行人工标注，人工标注方式包括更改文件名格式。

进一步地，所述步骤s1包括：由k-means控制机预先随机生成参数，所述参数包括聚类个数m和传入的微震信号序列x1-xn。

进一步地，所述步骤s2包括：

s21、通过dtw向量器对k-means控制机的所有微震信号序列x1-xn进行向量化；

s22、根据聚类个数确定初始聚类中心，并将初始聚类中心进行向量化；

s23、k-means控制机将每一个样本微震信号xi向量与初始聚类中心ti向量传给dtw剥离器；所述dtw剥离器将所述xi向量和ti向量构成一个距离矩阵yij，其中，距离矩阵yij中每一个元素的欧式距离d的计算公式如下：

d＝(xi,i-ti,j)²

式中，xii为xi向量的第i个点；tij为ti向量的第j个点。

进一步地，所述步骤s3包括：

s31、dtw剥离器将距离矩阵yij传给dtw距离计算器,dtw距离计算器在满足最佳匹配路径的前提下，计算累积距离y(i,j)，所述累积距离y(i,j)的计算公式如下：

y(i,j)＝d(xii,tij)+min{y(i-1,j-1),y(i-1,j),y(i,j-1)}

s32、k-means控制机得到所有样本微震信号与聚类中心之间的欧氏距离d；按照距离大小进行排序得到第一次聚类效果；

s33、k-means控制机将聚类后的第i类标签数据xi(xi1,xi2,xi3,xi4,…,xim),并令其中一个序列xim作为平均序列传输给dba更新器；

s34、通过与xim比较，使其他标签数据更新为新的序列，并将新的序列求平均作为新的聚类中心，再在dba更新器中循环上述操作，得到一个相对稳定的新的聚类中心ti，并将新的聚类中心ti返回给k-means控制机。

进一步地，所述最佳匹配路径包括：

a、从距离矩阵yij的起始点y(1,1)匹配到终止点y(n,m)；

b、距离矩阵yij中匹配后的前后两点一定是相邻关系；

c、在匹配过程中，假如点y(i,j)匹配下一个点，则下一个匹配点为y(i+1,j)、y(i,j+1)或y(i+1,j+1)中的一个，所述下一个匹配点选择方法为：分别计算y(i+1,j)、y(i,j+1)、y(i+1,j+1)与y(i,j)之间的欧氏距离d，选择欧氏距离d最小的那个点。

进一步地，所述步骤s4中，比较初始聚类中心与更新后的聚类中心是否相等包括：

k-means控制机将每一次dba更新器所得到的聚类中心都传递给聚类生成器进行存储；

聚类生成器得到下一个聚类中心ti与上一个聚类中心ti-1比较,判断两个聚类中心是否完全一致。

本发明提供的基于改进k-means的微震信号聚类装置，其包括：k-means控制机、dtw向量器、dtw剥离器、dtw距离计算器、dba更新器和聚类生成器，所述k-means控制机分别与所述dtw向量器、dtw剥离器、dtw距离计算器、dba更新器和聚类生成器相连接；其中，

dtw向量器，用于将从k-means控制机中输入的样本数据与聚类中心序列转换成向量化；

dtw剥离器，用于将每一个样本微震信号xi向量与初始聚类中心ti向量传构成一个距离矩阵yij；

dtw距离计算器，用于在满足最佳匹配路径的前提下，计算距离矩阵yij的累积距离y(i,j)；

dba更新器，用于根据聚类得到的结果更新新的聚类中心；

聚类生成器，用于输出聚类结果。

k-means控制机，分别与所述dtw向量器、dtw剥离器、dtw距离计算器、dba更新器和聚类生成器相连接。

本发明的有益效果在于：通过考虑微震信号具有时间序列的特点，采用dtw距离作为样本间的相似性度量，并改进k-means算法的聚类中心更新策略，采用与dtw特性相结合的dba质心更新方法，迭代优化聚类中心的更新过程，从而使其更适用于微震信号聚类分析，最大程度保留了波形的特征，其处理效率高、可靠性好。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明的基于改进k-means的微震信号聚类方法的聚类结果图；

图2是基于改进k-means的微震信号聚类方法的流程图；

图3是本发明基于改进k-means的微震信号聚类装置的结构示意图；

图4是本发明的一组不同特征的微震信号示例波形图；

图5是本发明实施例一的基于改进k-means的微震信号聚类装置的结构示意图；

图6是本发明实施例一的基于改进k-means的微震信号聚类方法的流程图；

图7是本发明实施例一的dtw剥离化向量示意图。

具体实施方式

下面结合附图对本发明的优选实施例进行详细阐述，以使本发明的优点和特征能更易于被本领域技术人员理解，从而对本发明的保护范围做出更为清楚明确的界定。

矿山微震是由矿山岩石破裂或流体扰动所引起的内部能量积累的震动。参阅图1所示，第一类为聚类装置产生的微震事件，第二类为在矿区生产过程中有机器、人为所产生的随机噪声。

参阅图2所示，本发明的基于改进k-means的微震信号聚类方法，其包括以下步骤：

s1、将微震信号传入到k-means控制机中，生成样本数据和初始聚类中心；

s2、计算样本数据与初始聚类中心之间的dtw距离；

s3、比较dtw距离大小，标注标签后进行第一次聚类，将第一次聚类得到的结果传给dba更新器中，得到新的聚类中心；

s4、比较初始聚类中心与更新后的聚类中心是否相等，若是，则执行步骤s5，若否，则执行步骤s6；

s5、输出聚类结果；

s6、返回步骤s3-s4，进行下一次聚类。

本发明通过考虑微震信号具有时间序列的特点，采用dtw距离作为样本间的相似性度量，并改进k-means算法的聚类中心更新策略，采用与dtw特性相结合的dba质心更新方法，迭代优化聚类中心的更新过程，从而使其更适用于微震信号聚类分析，最大程度保留了波形的特征，其处理效率高、可靠性好。

具体地，本发明的步骤s1之前包括步骤：

s0、对微震信号进行人工标注，人工标注方式包括更改文件名格式。这样可以更好地确认聚类结果。例如，以聚类结果来命名文件名的起始名，该命名方式只是为了最后得到聚类结果进行比较，而对k-means控制机进行聚类并无其他影响。

本发明的步骤s1包括：利用k-means的原理，由k-means控制机预先随机生成参数，参数包括聚类个数m和传入的微震信号序列x1-xn(x1,x2,x3,…,xn)。

本发明的步骤s2包括：

s21、通过dtw向量器对k-means控制机的所有微震信号序列x1-xn进行向量化；

s22、根据聚类个数确定初始聚类中心，并将初始聚类中心进行向量化；具体地，根据聚类个数m，采用人工指定的方法令m个样本微震信号作为初始聚类中心，这样可以提高微震信号的聚类平均准确率，初始聚类中心也是通过dtw向量器进行向量化。

s23、k-means控制机将每一个样本微震信号xi向量与初始聚类中心ti向量传给dtw剥离器；dtw剥离器将xi向量和ti向量构成一个距离矩阵yij，其中，距离矩阵yij中的每一个元素的欧式距离d的计算公式如下：

d＝(xi,i-ti,j)²

式中，xii为xi向量的第i个点；tij为ti向量的第j个点。

本发明的步骤s3包括：

s31、dtw剥离器将距离矩阵yij传给dtw距离计算器,dtw距离计算器在满足最佳匹配路径的前提下，计算累积距离y(i,j)，累积距离y(i,j)的计算公式如下：

y(i,j)＝d(xii,tij)+min{y(i-1,j-1),y(i-1,j),y(i,j-1)}

s32、k-means控制机得到所有样本微震信号与聚类中心之间的欧氏距离d；按照距离大小进行排序得到第一次聚类效果；

s33、k-means控制机将聚类后的第i类标签数据xi(xi1,xi2,xi3,xi4,…,xim),并令其中一个序列xim作为平均序列传输给dba更新器；

s34、通过与xim比较，使其他标签数据更新为新的序列，并将新的序列求平均作为新的聚类中心，再在dba更新器中循环上述操作，由于该循环过程是一个收敛的过程，所以最后得到一个相对稳定的新的聚类中心ti，并将新的聚类中心ti返回给k-means控制机。

本发明中的最佳匹配路径包括：

a、边界性准则：假设要匹配的两个向量的长度分别为(n,m)，需要从距离矩阵yij的起始点y(1,1)匹配到终止点y(n,m)；

b、连续性准则：距离矩阵yij中匹配后的前后两点一定是相邻关系(该相邻关系包括对角线关系)；使得整条匹配的最佳路径是一条连续的路径。

c、单调性原则：在匹配过程中，假如点y(i,j)匹配下一个点，则下一个匹配点为y(i+1,j)、y(i,j+1)或y(i+1,j+1)中的一个，下一个匹配点选择方法为：分别计算y(i+1,j)、y(i,j+1)、y(i+1,j+1)与y(i,j)之间的欧氏距离d，选择欧氏距离d最小的那个点。

本发明的步骤s4中，比较初始聚类中心与更新后的聚类中心是否相等包括：

k-means控制机将每一次dba更新器所得到的聚类中心都传递给聚类生成器进行存储；

在聚类生成器得到下一个聚类中心ti(ti,1,ti,2,…,ti,m)的时候与上一个聚类中心ti-1(ti-1,1,ti-1,2,…,ti-1,m)比较两个聚类中心是否完全一致，若两个聚类中心完全一致则聚类结束，聚类生成器生成聚类结果，并且将结果传输给k-means控制机。

参阅图3所示，本发明的基于改进k-means的微震信号聚类装置，其包括：k-means控制机、dtw向量器、dtw剥离器、dtw距离计算器、dba更新器和聚类生成器，k-means控制机分别与dtw向量器、dtw剥离器、dtw距离计算器、dba更新器和聚类生成器相连接；其中，

dtw向量器，用于将从k-means控制机中输入的样本数据与聚类中心序列转换成向量化；

dtw剥离器，用于将每一个样本微震信号xi向量与初始聚类中心ti向量传构成一个距离矩阵yij；

dtw距离计算器，用于在满足最佳匹配路径的前提下，计算距离矩阵yij的累积距离y(i,j)；

dba更新器，用于根据聚类得到的结果更新新的聚类中心；

聚类生成器，用于输出聚类结果。

k-means控制机，分别与dtw向量器、dtw剥离器、dtw距离计算器、dba更新器和聚类生成器相连接。

实施例一：

参阅图4-6所示，首先在矿区随机抽取200个微震信号，每一个微震信号的波形图类似于图4所示。每一个微震信号波形图实际上就是几千个采样点组合而成的一组时间序列。

分别将每个微震信号都通过dtw向量器转换为向量的形式存储起来，其存储形式为x1,x2,x3,…,xn。dtw向量器会将这些向量传回给k-means控制机存储。根据聚类个数k，指定k个样本微震信号向量作为初次聚类中心向量t1,t2,…,tk。

接着，分别将每个微震信号xi与每一个聚类中心tk通过dtw剥离器算出dtw距离矩阵，如图7所示，以x1,x2,x3,…,xn与聚类中心t1为例：分别得到距离矩阵y11,y21,…,yn1。以矩阵y11为例，矩阵中的每一个点y11(i,j)均是微震信号x1的第i个值与聚类中心t1中的第j个值之间的欧式距离d，这里假如微震信号x1为{1,2,1,3,4},聚类中心t1为{4,2,6,2,1}，根据dtw计算距离的原理，这里所得到的y11矩阵形式为：

在得到所有的距离矩阵y之后，都会将矩阵传入dtw距离计算器，通过匹配最佳路径的方法，算出相应的样本数据与聚类中心之间的dtw距离。

以y11矩阵为例：从y11的左上角第一个元素开始进行匹配，在临近元素中{0,1,4}选择一个最小的进行匹配，此时累加距离v11＝9+0＝9。接下来，从该位置开始继续匹配，在临近元素中{1,25,16}中进行匹配，则此时累加距离v11＝9+0+1＝10。接下来从该位置开始继续匹配，在临近元素中{1,25,9}中进行匹配，则此时累加距离v11＝9+0+1+1＝11。接下来从该位置开始继续匹配，在临近元素中{4,4,9}中进行匹配，则此时累加距离v11＝9+0+1+1+4＝15。由于这里遵循着最佳匹配路径，所以选择最近的路径，即对象线的位置{4}。一直匹配到最后得到的累加距离为v11＝9+0+1+1+4+4+9＝28。根据dtw原理，匹配出最佳路径，以y’11所示：

用以上的方法所得到的所有的dtw距离传入k-means控制机，按就近原则进行第一次聚类分类。在dtw所得到的最佳匹配关系的基础上形成新的时间序列，在此基础上求得簇内平均序列，以平均序列为新的类内中心，簇内样本重复上述操作，通过迭代来更新聚类中心。

本发明的dba更新的思想如下：

(1)根据dtw的匹配最佳路径的原则，所得到的a，b，c序列为dtw距离计算器得到的最佳匹配序列，若a为初始指定平均序列，b和c为簇内样本序列。

首先根据最佳匹配思想进行分别将b和c与a进行匹配。

(2)得到最佳的匹配关系后，根据一对多和多对一的思想分别对其进行合并和分裂，由于a为平均序列，所以a不发生变化，仅仅改变b和c的序列。

这里以a与c匹配距离进行阐述，在匹配前a中<2>对应c中的<1,1>，所以根据一对多的原则，需要将c中的<1,1>进行合并元素得到<1>，而对于a中的<0,0>所对应的c中<3>，根据一对多的原则将c中的<3>分裂得到<3,3>。

(3)这里更新了b和c序列以后，接下来进行更新平均序列(这里的平均序列为聚类中心)的操作，即对所有序列对应元素求平均值得到新的平均序列arg＝<1/3,4/3,21/3,11/9,1>。接下来重复1和2中的操作分别将a，b和c与arg进行匹配更新平均序列，直到聚类中心收敛不再变化。

a＝<0,2,10,0,0>

b＝<1,1,9,2/3,0>

c＝<0,1,2,3,3>

arg＝<1/3,4/3,21/3,11/9,1>

(4)此时，第一次聚类结束，但是聚类效果还未达到要求，所以得到新的聚类中心后，再将新的聚类中心放入k-means控制机中重复循环，最开始的操作，直到dba更新器中所得到的聚类中心不发生变化时，标志着聚类结束，此时由聚类生成器输出聚类结果。

以上，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何不经过创造性劳动想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书所限定的保护范围为准。