一种基于蚁群算法改进的LDTW序列相似度量方法与流程

本发明属于机械故障诊断和分类领域，其实质为一种基于蚁群算法改进的ldtw序列相似度量方法。

背景技术：

相似度量通常作为数据挖掘的前序步骤，与挖掘算法结合起来完成不同的挖掘任务。挖掘算法在训练过程中，通常都需要计算序列间的相似度，发现它们之间的相似规律或特征，许多研究表明，经过精心设计的相似性度量可以提升算法的性能。动态时间弯曲(dynamictimewarping，dtw)作为如今很常用的序列间的相似度量方法,可以在时间或空间发生扭曲变形的情况下进行相似度量。目前，dtw的应用已经从语音识别，发展到其他诸多领域中，如生物信息学，在线签名验证和手势识别等。

如今，随着数据量的增多，dtw计算所花费的时间也随之增加；同时，dtw还存在病理对齐的问题，即一个序列上的一个点可能会映射到另一个序列的多个点上，形成多个奇异点，从而导致dtw的精度也有所下降。

由于dtw在时间和空间上的复杂性受到了限制，众多学者对其性能和效率上进行了改进；其中包括很多经典的改进算法，大致将这些变体分为两大类：

第一类是取代了对dtw特征的考虑，研究者觉得数据的趋势更为重要。例如，微分dtw(ddtw)通过考虑序列的一阶导数来获得形状的信息，从而减少奇异点的发生，但是该方法没有考虑算法的时间和空间的复杂性。定性dtw(qdtw)考虑了dtw时间复杂度高的问题，将时间序列转换为定性时间序列，从而来减少时间序列的长度，再计算定性时间序列之间的dtw，但是实验中发现预测精度略有下降。还有很多研究是针对其他数据特征的考虑，例如形状上下文和兴趣点等。本发明就不再赘述了。

第二类基本上是对弯曲路径设置约束的。其中，分割dtw(sdtw)证明了用分段线性段支持时间弯曲查询的可能性，但是仍存在病理对齐问题。加权dtw(wdtw)对参考点与测试点之间相位差较高的点进行处罚，以防止异常值造成的最小距离失真，实验结果表明提高了分类和聚类精度，但是计算所花费的时间和空间是很大的。快速搜索dtw(fdtw)考虑了dtw计算时间复杂度高，采用有效地删除大量的搜索候选项，从而直接降低搜索成本。但是其中采用了近似的方法，用几个片段表示一个序列的粗略版本来估计时间弯曲距离，很可能错过了正确的对齐。量化dtw采用了码字的概念，在做语音识别的时候，不再存储同一个单词的多个模板，而是为每个单词存储一个参考模型，将未知语音与参考模型的中心体进行比较，提高了识别的速度，但是识别率不可避免的有所下降。柔性dtw(fdtw)给予额外的分数来奖励较长的一对一的片段，避免了一个序列中的相邻点与另一个序列中相隔较远的点对齐。由于额外的增加项导致了其计算时间的增加。有限弯曲长度下的动态时间弯曲距离(ldtw)是最近提出的改进dtw的方法，属于第二类对弯曲路径提出一种新的约束，但该方法无法在保证精确性的前提下有着较低的时间复杂度。

这些方法要么取代了对dtw特征的考虑要么对弯曲路径设置约束，但都无法做到在保证精确性的前提下尽可能的降低时间复杂度。

技术实现要素：

有鉴于此，本发明提出了一种基于蚁群算法(antcolonyoptimization,aco)改进的ldtw的序列相似度量方法，即aco_ldtw；本发明不仅能有效地抑制病理对齐，同时允许更多的灵活性。如图1所示，其具体包括以下步骤：

s1、将给定长度的两个序列转换为距离矩阵；

s2、根据ldtw距离的最大步长确定限制距离，将距离矩阵根据限制距离映射到栅格法的0-1矩阵中，根据限制距离确定的搜索范围，按照禁忌表对蚂蚁在栅格中的搜索范围进行限制；

s3、蚂蚁限制后的栅格区域移动，根据动态调整状态转移概率移动到下一访问点；从而搜索得到路径；

s4、根据蚂蚁搜索得到的栅格值之和，确定适应度函数；在适应度函数最小时，即为最佳的序列距离长度；

s5、根据最佳的序列距离长度，确定两个序列的相似度；

其中，ldtw距离表示为有限弯曲长度下的动态时间弯曲距离。

进一步的，在给定长度分别为i、j的两序列x＝{x1,x2,…,xi}和y＝{y1,y2,…,yj}中，称

为序列x和y的距离矩阵，其中元素dis(xi,yj)为序列x中第i个观测数据和序列y中的第j个观测数据之间的距离，表示为dist(xi,yj)＝(xi-yj)²。根据对齐矩阵，可将序列x与y间的dtw距离做如下定义：

进一步的，基于栅格地图的蚁群算法采用的0-1矩阵模拟环境(“0”代表无障碍物的栅格，“1”代表有障碍物的栅格)，通过禁忌表限制蚂蚁已访问过的位置，在避开障碍物的情况下，搜索从起点到终点的最优路径。

具体的，所述步骤s2包括根据限制距离，将距离矩阵中每个元素映射到栅格法的0-1矩阵中；即限制距离将蚂蚁在栅格中的搜索范围限制为六边形，所述六边形的两个对角为直角，直角边的距离均为限制距离l，非直角边上的栅格为“1”，即该栅格的矩阵元素值为“1”；其余栅格均为“0”；蚂蚁避开为“1”的栅格，且避开禁忌表上已访问过的栅格，搜索从起点到终点的最优路径；

其中，限制距离表示为l＝ceil[(nstep_max-ndim)×ε]，nstep_max为ldtw距离的最大步长，ndim是序列的维数，ceil[]表示向正无穷方向取整，ε表示限制因子，ε∈[0.1,0.4]；“0”代表无障碍物的栅格，“1”代表有障碍物的栅格。如图2所示，该图的可搜索范围被黑色障碍限制为六边形。

同时，ldtw引入步长限制了弯曲路径范围，避免了病理性对齐问题。通过公式可以看出，ldtw是通过从点(i,j)出发进行递归计算到(1,1)，且步数不能超过所限制的最大步长。本发明将对齐矩阵映射成栅格地图中的0-1矩阵，再利用本发明提出的算法搜索最优路径。当数据维数过大时，ldtw的计算效率会显著降低，本发明为提升计算效率，结合ldtw的步数范围，将蚂蚁在地图中的搜索范围进行限制。

系数ε控制着搜索范围的范围，该系数的增减会导致搜索范围的增减。为了讨论该系数的取值范围，选用了本发明实施例1中第一组维数为60×60的syntheticcontrol数据，其最大步数为72，分别将ε设为不同的值，用本发明算法进行搜索，在计算迭代100次后，得到的效果如图3-6所示。

图2中，可搜索范围被黑色栅格限制在中间的狭窄区域内，图3中，曲线路线表示一个搜索结果。其中，ε分别取0.4,0.3，0.2，0.1。可以看出，当ε取值较大时，搜索范围过大，算法不能有效的收敛到最优解，且得到结果中，步长已经超过了最大步数；而ε取值过小时，搜索范围变小，导致遗漏了部分优质解。因此本发明将ε取值为0.2，既符合ldtw的步数限制，又通过缩小搜索范围提高了蚁群算法计算效率。

进一步的，考虑到ldtw的三个约束规则，本发明将蚂蚁的每一次行走区域方向进行了限制。如图7所示，蚂蚁在选择下一个访问点时最多只能在三个栅格中选择，传统蚁群算法常采用如下的状态转移概率来选择下一个访问点；也即所述步骤s3包括蚂蚁在限制后的限制区域进行移动，通过生成一个均匀随机数rand()和对比参数p，如果rand()>p，蚂蚁在当前栅格用原始状态转移概率公式从可访问栅格中选择下一个访问栅格，否则，在可访问的栅格中选择较小的值为下一个访问栅格，其中，可访问的栅格为当前栅格相邻的三个栅格，分别为其右上方栅格、右方栅格以及上方栅格。

进一步的，所述原始状态转移概率公式为：

其中，表示蚂蚁k在第t次迭代中走(i,j)栅格时的概率；τij(t)表示第t次迭代时(i,j)栅格的信息素，τis(t)表示第t次迭代时(i,s)栅格的信息素；ηij(t)表示在第t次迭代时栅格(i,j)的启发式因子，ηis(t)表示第t次迭代时栅格(i,s)的启发式因子，在本发明中，启发式因子为两个序列中相应点之间距离的倒数；常数α表示信息启发因子，常数β表示期望启发因子；allowedk表示蚂蚁的下一个访问栅格。

优选的，所述第t次迭代时的信息素是经过如下处理：

当t＝0时，对信息素进行初始化，将对角线上的每个栅格的信息素初始化为常数c，在搜索范围内的其余栅格的信息素则在[0,c]中随机取值；

当t≠0时，当第t次迭代后的新路径优于当前最优路径时，则采用混合强化公式对(i,j)栅格的信息素更新；否则，根据最大-最小蚂蚁系统中的更新公式对(i,j)栅格的信息素进行更新；

其中，混合强化公式表示为：

τij表示栅格(i,j)的信息素；表示走到栅格(i,j)时的信息素总数，表示第(i,j)栅格信息素的增量，其计算公式如下：

ft为适应度评价函数，sibest为迭代后的最优解集，即最优路径。

进一步的，在进行第t次迭代更新时，若信息素出现停滞，则根据平滑公式对信息素进行更新，再利用混合强化公式对信息素更新；平滑公式表示为：τij→ρ×(τmax-τij)；τ1是指停滞不前的路径上的信息素，τ2表示由平滑机制处理过后的信息素，τmax表示最大信息素值；ρ表示信息素的蒸发系数，取ρ＝0.1。

进一步的，由于两者算法(蚂蚁算法和ldtw距离算法)在应用背景上有一定的差异，基于栅格地图的aco中适应度函数通常是寻找起点到终点的最短路径，但ldtw是找到最小累计成本矩阵中的最佳弯曲距离。因此本发明在栅格搜索路径时，将路径的长度转化为对应弯曲距离，适应度函数(弯曲距离)，也即是步骤s4的适应度函数为：

其中，sk表示蚂蚁k从起始栅格走到目标栅格构成的一个解；dis(q)表示的是蚂蚁k所走在路径中第q个格子中的值，该值对应为序列x和序列y中相应点之间的距离；1<n≤lub+1；lub是ldtw中弯曲路径的上界；n表示该路径对应的总格子数。

本发明将原始的蚁群算法中的状态转移概率改进为动态调整状态转移概率。通过生成一个0-1均匀随机数rand()和参数p，如果rand()>p，则用原始状态转移概率公式选择下一个访问点，否则，在可以访问的栅格中选择较小的值为下一个访问节点。参数p随着时间的推进而递增，本发明主要采用蝙蝠算法中的递增算子：

p^t+1＝p⁰[1-e^-γt]

式中，p^t+1表示在t+1次迭代时的对比参数；γ为常数，取值为[0.85,0.95]；p⁰表示初始对比参数，为[0,1]间的任意数，需要注意的是，随着计算时间的流逝，参数p会越来越接近p⁰。通过该操作保证了前期搜索具有随机性，而在后期加快了收敛速度。

进一步的，由于蚁群算法本身也存在一定缺陷，导致收敛速度慢，已陷入局部最优等问题。这些缺陷主要体现在信息素的相关操作中。首先在信息素的初始化阶段，传统的蚁群算法，如最大-最小蚂蚁系统(mmas)将图中每个位置的信息素含量初始化为常数c，导致前期搜索速率下降，因此本发明为了让算法前期尽量在搜索范围的对角线上搜索，将该对角线上的信息素初始化为常数c，在搜索范围内的其余栅格的信息素则在[0,c]中随机取值。

具体的，在算法运行时蚁群更倾向于选择信息素含量高的区域，易忽略信息素含量较低的区域。尤其是原始蚁群算法的更新操作会使得不同区域的信息素含量差异越来越大，导致算法陷入局部最优，甚至出现停滞现象。为解决这一问题，本发明采用jiaxuning提出的一种最优路径更新操作，该方法为了减少算法陷入局部最优，引入了一种混合强化信息素更新机制和信息素平滑机制来解决原始算法的停滞现象，具体操作如下：首先，对当前信息素矩阵中停滞路径上的信息素值进行如下处理：

τij→ρ×(τmax-τij)

τij是指停滞不前的路径上的信息素，ρ表示信息素的蒸发系数，取ρ＝0.1。若算法在计算中期，出现停滞现象，这利用上式对该路径进行平滑处理。当迭代中两个解集的相识度超过了95％就认定处于停滞状态，即在第t次和第t-1次迭代中，两者解集的相似度超过95％就认定信息素停滞。

然后，引入混合强化信息素更新机制，对信息素进行更新：

τij表示栅格(i,j)的信息素；表示走到栅格(i,j)时的信息素总数；表示第(i,j)栅格信息素的增量，它增强了新探索的路径。

在迭代过程中，当新路径优于当前最优路径时，在混合强化方程中加入否则，根据最大最小蚂蚁系统(mmas)中的信息更新方程进行信息素更新。计算公式如下：

ft为适应度评价函数，sibest为迭代后的最优解集，即最优路径。当算法接近停滞时，使用新的信息素平滑机制重新初始化信息素矩阵以平衡信息素浓度。以平衡和路径的信息素浓度差，这有助于蚂蚁在接近停滞的时期继续搜索。

本发明提出了一种新的计算序列相似度量的方法，称为aco_ldtw。通过对适应度函数，状态转移概率，信息素初始化和更新机制的改进，同时将蚂蚁搜索范围进行了缩减，既保证了弯曲路径长短不会超过ldtw限制的弯曲路径长度，又降低计算的时间成本。对机械故障诊断和分类领域有实质性的帮助。

本发明的有益效果：

本发明提出了一种新的计算序列相似度量的方法。通过对适应度函数，状态转移概率，信息素初始化和更新机制的改进，同时将蚂蚁搜索范围进行了缩减，既保证了弯曲路径长短不会超过ldtw限制的弯曲路径长度，又降低计算的时间成本。对机械故障诊断和分类领域有实质性的帮助。

附图说明

图1为本发明的流程图；

图2为本发明采用的对齐矩阵映射的栅格地图；

图3为本发明经过计算迭代后ε对搜索范围的影响示意图a；

图4为本发明经过计算迭代后ε对搜索范围的影响示意图b；

图5为本发明经过计算迭代后ε对搜索范围的影响示意图c；

图6为本发明经过计算迭代后ε对搜索范围的影响示意图d；

图7为本发明基于栅格地图蚁群算法的最优路径；

图8为未采用本发明的原始rbf核函数的结果；

图9为采本发明后的核函数结果；

图10为本发明与欧式距离方法的1nn应用错误率比较；

图11为本发明与dtw方法的1nn应用错误率比较；

图12为本发明与dtw(bestsakoe-chibaband)方法的1nn应用错误率比较；

图13为本发明与dtw(bestltakuraparallelogram)方法的1nn应用错误率比较；

图14为本发明与ldtw方法的1nn应用错误率比较；

图15为本发明与dtw、ldtw方法在六个不同数据集上进行1nn分类的运行时间。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。

本发明提供了一种基于修正后的蚁群算法的序列相似度量方法，如图1所示，它包括如下步骤：

s1、将给定长度的两个序列转换为距离矩阵；

s2、根据ldtw距离的最大步长确定限制距离，将距离矩阵根据限制距离映射到栅格法的0-1矩阵中，根据限制距离确定的搜索范围，按照禁忌表对蚂蚁在栅格中的搜索范围进行限制；

s3、蚂蚁限制后的栅格区域移动，根据动态调整状态转移概率移动到下一访问点；从而搜索得到路径；

s4、根据蚂蚁搜索得到的栅格值之和，确定适应度函数；在适应度函数最小时，即为最佳的序列距离长度；

s5、根据最佳的序列距离长度，确定两个序列的相似度；

其中，ldtw距离表示为有限弯曲长度下的动态时间弯曲距离。可以理解的是，本发明中步骤s5根据序列距离长度确定两个序列的相似度，是本领域技术人员所熟知的技术，本发明对此不做具体的限定。

实施例1

在凯斯西储大学的轴承故障数据中，数据是利用电火花加工切槽来模拟轴承不同部位的单点故障，主要是通过设置不同的切槽宽度：0.007、0.014、0.021(1英寸＝2.54厘米)，来模拟轴承轻度故障、中度故障以及严重故障。本次实验选用主轴转速1772rpm、载荷为1hp、采样频率为12000hz的驱动端振动轴承的6种故障振动信号进行分析，包括：外圈轻微故障、内圈轻微故障、滚动体轻微故障、内圈中度故障、内圈严重故障以及滚动体严重故障。具体的步骤可如下：

step1:在每种故障状态下测取50组样本，每个样本采样点数为2048点，将其中20组样本作为训练样本，其余30组样本作为测试样本；

step2：从每个样本的振动信号中提取24个时频域统计特征；

step3:采用支持向量机(svm)对以上数据进行故障诊断。首先采用高斯核函数(rbf)作为svm中的核函数，其表达式如下：

step4:上式中用||xi-xj||²计算样本间的相似距离，为验证本发明算法的有效性，将||xi-xj||²替换为aco_ldtw(xi,xj)计算所得的相似距离；aco_ldtw为本发明提出的基于蚁群算法改进的ldtw序列相似度量方法，为了方便描述，以下都简称为aco_ldtw。

step5：由于支持向量机的核函数和惩罚函数对其分类性能影响较大，本发明选用蝙蝠算法对支持向量机参数进行调优，分别选用原始的高斯核函数(rbf，径向基函数)和替换后的核函数(aco_ldtw方法的核函数)进行调参，以得到两者最优的分类准确率。两种核函数的分类准确率收敛曲线如图8，9所示：

由图8，9可以看出，采用aco_ldtw(xi,xj)改进的核函数能在一定程度上提升准确率，平均准确率达到了99.1667％，optimalfitness表示最佳适应度值，averagefitness表示平均适应度值，accuracy表示准确率；因此本发明在实际工程问题中具有一定的应用能力。

实施例2

为了有效解决ldtw的计算时间过长的问题，因此本节中，在文献1nn分类实验的基础上评估aco_ldtw作为距离度量的性能。

数据集是采用ucr时间序列分类存档的二十二个数据集上进行的。它主要包括从各种领域收集的实时序列数据集和一些合成的数据集。表1总结了实验中使用的所有数据集的基本信息，如类的数目、训练集的大小、测试集的大小和时间序列长度。类的数量从2(如gun-point、wafer等)到50(50words)。时间序列的长度范围从60(syntheticcontrol)到637(lightning-2)。

其中，aco_ldtw的相关参数设置如下：α表示状态转移概率公式的信息素部分的权重，通常设置为2；β表示启发式因子的加权值，通常设置为8。ρ为信息素挥发系数，通常在[0，1)内取值，根据文献中的参数实验，本发明将ρ设置为0.02。

一个最近邻分类(1nn简称)被强烈推荐来证明距离测度的性能，因为1nn没有参数，因此精度仅仅依赖于距离度量。ldtw的1nn错误率与最佳的翘曲路径长度上限文献都已被给出。表1给出了ucr分类文档中22个数据集的基本信息：

表1ucr分类文档中22个数据集的基本信息

本发明将aco_ldtw与ldtw，还有四个最广泛使用的距离度量方法：euclideandistance，naivedtw，dtwunderbestsakoe-chibaband和dtwunderbestitakuraparallelogram进行对比。从表1中，可以观察到aco_ldtw在22个数据集中有16个数据集的错误率最低，除了gun-point，face(all)和lightning-7，其他数据集都与ldtw在1nn错误率上保持一致。并且比ldtw花费的计算时间更短，这些结果证明了所提出的aco_ldtw的实用性。图10-14示出了aco_ldtw与其它五个方法之间的1nn错误率的两两比较。在这样的设置中，其中一条折线是由aco_ldtw在不同数据集上的1nn错误率形成的，另一条代表其他对比方法。其中，离圈中心越近的点表示其错误率越小，且优于在该数据集上的其他方法，反之亦然。在图10-14中，aco_ldtw对应的大部分点距离圈中心点越近，说明其效果就越好。在每个图上描绘了aco_ldtw更靠近圈中心和其他对比方法更靠近圈中心的数量，其余为两者距离圈中心相等的个数。例如，在图10上，aco_ldtw在22个数据集中的有18个数据集优于欧氏距离euclideandistance，在4个数据集上相等，并且没有在任何数据集上劣于euclideandistance。同理如图11，aco_ldtw在22个数据集中的有15个数据集优于dtw，在6个数据集上相等。这些图表直观地证明aco_ldtw在各种数据集上的优势。其中，如表2所示，本发明与其他方法的错误率进行比较；

表21nn错误率比较

接下来，对aco_ldtw的时间复杂度进行分析和运行时间与其他方法进行比较。lub是aco_ldtw中弯曲路径的上界。aco_ldtw的迭代次数为k，蚂蚁数量为m的时间复杂度o(k·m·lub)。

为了比较不同方法的运行时间，将五个不同数据集在进行1nn分类实验基础上记录了相应的计算时间，图15显示了平均运行时间。本发明可以看到，aco_ldtw比ldtw的计算时间更短。因为采用动态调整状态转移概率保证了前期搜索具有随机性，而在后期加快了收敛速度。并且蚂蚁不用在整个栅格地图中寻找最优路径，而是通过限制了搜索范围，范围越小，计算时间越短，同时也节约了计算时间。

在数据挖掘中，本发明提供的基于蚁群算法改进的ldtw的序列相似度量方法对机械故障诊断和1nn分类能够提供更准确高效的结果，可以通过本发明研究数据挖掘任务提供帮助。

本领域普通技术人员可以理解上述实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：rom、ram、磁盘或光盘等。

以上所举实施例，对本发明的目的、技术方案和优点进行了进一步的详细说明，所应理解的是，以上所举实施例仅为本发明的优选实施方式而已，并不用以限制本发明，凡在本发明的精神和原则之内对本发明所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。