CN111189638B - 基于hmm和qpso优化算法的轴承故障程度辨识方法 - Google Patents

基于hmm和qpso优化算法的轴承故障程度辨识方法 Download PDF

Info

Publication number
CN111189638B
CN111189638B CN201911342148.2A CN201911342148A CN111189638B CN 111189638 B CN111189638 B CN 111189638B CN 201911342148 A CN201911342148 A CN 201911342148A CN 111189638 B CN111189638 B CN 111189638B
Authority
CN
China
Prior art keywords
formula
follows
hidden markov
iteration
algorithm
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201911342148.2A
Other languages
English (en)
Other versions
CN111189638A (zh
Inventor
党鹏飞
王明罡
杨铮鑫
鲍宁波
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shenyang University of Chemical Technology
Original Assignee
Shenyang University of Chemical Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shenyang University of Chemical Technology filed Critical Shenyang University of Chemical Technology
Priority to CN201911342148.2A priority Critical patent/CN111189638B/zh
Publication of CN111189638A publication Critical patent/CN111189638A/zh
Application granted granted Critical
Publication of CN111189638B publication Critical patent/CN111189638B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/04—Bearings
    • G01M13/045—Acoustic or vibration analysis

Landscapes

  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • General Physics & Mathematics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

基于HMM和QPSO优化算法的轴承故障程度辨识方法，涉及滚动轴承故障诊断方法，本发明通过对原始时态的信号作为输入，采用变分模态分解（VMD）的方法对信号进行分解，然后对每一个本征模态分量（IMF）分量中进行奇异值分解，得到奇异值矩阵，利用K‑means聚类方法，将奇异值矩阵分类并编码，最后利用量子粒子群算法（QPSO）进行对隐马尔科夫模型（HMM）的参数优化，得到训练好的HMM，再用同样的方式处理测试信号，然后用训练好的HMM对测试信号进行精确的故障程度辨识。本发明用参数作为隐马尔科夫模型的学习参数，从而对轴承故障程度进行分类，用以故障程度评估，具有重要的工程意义。

Description

基于HMM和QPSO优化算法的轴承故障程度辨识方法
技术领域
本发明涉及一种滚动轴承故障辨识方法，特别是涉及一种基于HMM和QPSO优化算法的轴承故障程度辨识方法。
背景技术
滚动轴承是旋转机械中的关键部件，其故障可能导致高成本停机，甚至造成整个机械的灾难性故障。为了保证机械的运行安全，降低维修成本，以故障程度评估技术为核心的状态维修越来越受到人们的重视。所以说对滚动轴承进行故障程度评估具有重要的工程意义。
目前常用的轴承故障程度辨识方法是对信号的时域或频域进行特征提取，然后通过分类器进行诊断，用隐马尔科夫模型进行故障程度辨识成为了一种很好地方法。但是这种方法具有局限性，通过隐马尔科夫模型的时间序列分类能力的确可以很好的将故障进行分类，但是用于隐马尔可夫模型训练的Baum-Welch算法需要大量的数据对其进行训练，而且可能会进入局部最优，无法得到最优化的隐马尔科夫模型参数，从而可能使分类的结果不准确。因此，本发明提供一种新的隐马尔科夫模型参数的优化方法，用量子粒子群算法（QPSO）面向全局搜寻隐马尔科夫模型的最优参数解，克服Baum-Welch算法的局部最优问题。
发明内容
本发明的目的在于提供一种基于HMM和QPSO优化算法的轴承故障程度辨识方法，通过变分模态分解（VMD）对轴承振动信号进行分解，然后使用奇异值分解（SVD）进行特征提取，并用K-means进行聚类，最后用量子粒子群算法（QPSO）面向全局搜寻隐马尔科夫模型的最优参数。用此参数作为隐马尔科夫模型的学习参数，从而对轴承故障程度进行分类。
本发明的目的是通过以下技术方案实现的：
基于HMM和QPSO优化算法的轴承故障程度辨识方法，所述方法包括如下步骤：(1)、采集滚动轴承振动信号使用振动数据采集仪以采样频率为12000Hz采集待检测滚动轴承在运行状态下的振动信号，记为X[m]，m＝1 ,2 ,…,M，M为总采样点数，并标记轴承状态，总状态数为y，将所有的所得的数据集合为，标签数据集为；(2)、数据的分帧操作，建立训练集取前120k个采样点，平均分成12组数据，每组数据10k个采样点，再把每组数据平均分成5帧，每帧2k个采样点；(3)、振动信号的分频处理利用变分模态分解VMD的可变尺度分解对训练集中的信号进行处理，变分模态分解的过程是将一个复杂的滚动轴承时域振动信号分解为K个IMF分量，K为任意正整数；
(4)、对本征模态分量进行奇异值分解（SVD）；利用奇异值分解（SVD），提取由变分模态分解（VMD）分解后的最优IMF的奇异值矩阵；
(5)、对奇异值向量进行聚类；
(6)、利用量子粒子群算法对隐马尔科夫模型的参数
Figure DEST_PATH_IMAGE002
进行估计；每一个粒子都是一个HMM参数，令其适应度函数为：
Figure DEST_PATH_IMAGE004
(13)；
(7)、诊断轴承；
(8)、诊断结束。
所述的基于HMM和QPSO优化算法的轴承故障程度辨识方法，所述变分模态分解的基本过程如下：
(3.1)、希尔伯特变换，对每一个IMF分量进行希尔伯特变换，得到单边频谱如下：
Figure RE-DEST_PATH_IMAGE005
（1）
式中
Figure RE-RE-DEST_PATH_IMAGE006
是第k个IMF分量，k为任意正整数，
Figure RE-569622DEST_PATH_IMAGE006
可看做谐波信号：
Figure RE-DEST_PATH_IMAGE007
，其中
Figure RE-RE-DEST_PATH_IMAGE008
是
Figure RE-237495DEST_PATH_IMAGE006
的瞬时幅值，且
Figure RE-DEST_PATH_IMAGE009
为
Figure RE-892598DEST_PATH_IMAGE006
的相位，且
Figure RE-RE-DEST_PATH_IMAGE010
，
Figure RE-DEST_PATH_IMAGE011
为冲激函数；
Figure RE-RE-DEST_PATH_IMAGE012
，t是时间；
Figure RE-DEST_PATH_IMAGE013
为卷积运算；
(3.2)、频率混合，利用上述得出的单边频谱信号混合一个估计的中心频率
Figure RE-RE-DEST_PATH_IMAGE014
，待混合之后，每一个
Figure RE-913906DEST_PATH_IMAGE006
的频谱将移动到相应的基频带上，即：
Figure RE-DEST_PATH_IMAGE015
（2）
(3.3)、带宽估计，计算公式（2）所得到的解调信号的梯度二范数平方：
Figure RE-RE-DEST_PATH_IMAGE016
（3）
其中
Figure RE-DEST_PATH_IMAGE017
为函数对时间t的偏导数；
(3.4)、建立最优化模型，通过以上计算，可以得到的变分约束模型如下：
Figure RE-RE-DEST_PATH_IMAGE018
（4）
式中，K为IMF分量个数；
Figure RE-DEST_PATH_IMAGE019
是模态的集合，
Figure RE-RE-DEST_PATH_IMAGE020
是
Figure RE-253181DEST_PATH_IMAGE019
的中心频率；
Figure RE-DEST_PATH_IMAGE021
是K个IMF
Figure RE-826376DEST_PATH_IMAGE019
的加和；
(3.5)、转化为非约束问题，引入二次惩罚因子
Figure RE-RE-DEST_PATH_IMAGE022
和拉格朗日乘子
Figure RE-DEST_PATH_IMAGE023
构成表达式：
Figure RE-RE-DEST_PATH_IMAGE024
（5）
Figure RE-844142DEST_PATH_IMAGE019
，
Figure RE-43042DEST_PATH_IMAGE020
和
Figure RE-DEST_PATH_IMAGE025
初始值都为0，通过交替方向乘子法（ADMM）迭代得到（4）的最优解，也可以得到
Figure RE-587287DEST_PATH_IMAGE019
所对应的中心频率
Figure RE-456017DEST_PATH_IMAGE020
如下所示：
Figure RE-RE-DEST_PATH_IMAGE026
（6）
Figure RE-DEST_PATH_IMAGE027
（7）
式中：
Figure RE-RE-DEST_PATH_IMAGE028
为当前余项
Figure RE-DEST_PATH_IMAGE029
的维纳滤波，
Figure RE-RE-DEST_PATH_IMAGE030
为模态功率谱的重心；
Figure RE-DEST_PATH_IMAGE031
为
Figure RE-RE-DEST_PATH_IMAGE032
的傅里叶变换，n为任意正整数，
Figure RE-DEST_PATH_IMAGE033
为频率，i为任意正整数；
(3.6)、变分模态分解（VMD）对滚动轴承轴承振动信号分解过程如下：设置n=n+1，k=k+1，依照公式（6）（7）更新
Figure RE-587177DEST_PATH_IMAGE019
和
Figure RE-324189DEST_PATH_IMAGE020
，直到k=K是结束迭代；依照公式（8）更新
Figure RE-RE-DEST_PATH_IMAGE034
公式如下：
Figure RE-DEST_PATH_IMAGE035
（8）
式中
Figure RE-RE-DEST_PATH_IMAGE036
为噪声容限参数，通常令
Figure RE-DEST_PATH_IMAGE037
，重复上述操作直到满足公式（9），就可以得到由变分模态分解（VMD）处理后的最优化IMF；公式如下：
Figure RE-RE-DEST_PATH_IMAGE038
（9）
式中
Figure RE-DEST_PATH_IMAGE039
为误差，可以取任意正整数。
所述的基于HMM和QPSO优化算法的轴承故障程度辨识方法，所述奇异值分解公式如下：
Figure 100002_DEST_PATH_IMAGE080
（10）
其中
Figure DEST_PATH_IMAGE082
为左奇异向量矩阵，
Figure DEST_PATH_IMAGE084
为右奇异向量矩阵，
Figure DEST_PATH_IMAGE086
为奇异值向量矩阵，它是一个对角矩阵，其主对角线上的值为矩阵IMF的奇异值，其他位置上的元素值都为0；用此方法，分别求出每组数据的奇异值向量。
所述的基于HMM和QPSO优化算法的轴承故障程度辨识方法，所述对奇异值向量进行聚类，使用K-means聚类算法，将每组数据进行分类，其算法流程如下：（1）从数据集中随机选取k个样本作为初始的k个质心向量：
Figure DEST_PATH_IMAGE088
；
（2）对于
Figure DEST_PATH_IMAGE090
（a）将簇划分
Figure DEST_PATH_IMAGE092
初始化为
Figure DEST_PATH_IMAGE094
；
（b）对于
Figure DEST_PATH_IMAGE096
，计算样本
Figure DEST_PATH_IMAGE098
和各个质心向量
Figure DEST_PATH_IMAGE100
的距离，公式如下：
Figure DEST_PATH_IMAGE102
（11）
将
Figure DEST_PATH_IMAGE104
标记最小的
Figure DEST_PATH_IMAGE106
所对应的类别
Figure DEST_PATH_IMAGE108
，此时更新
Figure DEST_PATH_IMAGE110
；
（c）对于
Figure DEST_PATH_IMAGE112
，对
Figure DEST_PATH_IMAGE114
中所有的样本点重新计算新的质心，公式如下：
Figure DEST_PATH_IMAGE116
（12）
（d）如果所有的k个质心向量都没有发生变化则转到步骤（c）
（3）输出簇划分；由此得到了一组整数，作为隐马尔可夫模型的观测数据，用于后期的故障辨识。
所述的基于HMM和QPSO优化算法的轴承故障程度辨识方法，所述利用量子粒子群算法求隐马尔科夫模型的最优化参数步骤如下：(6.1)、初始化，种群数量10，最大迭代次数30，算法运行10轮；(6.2)、随机生成隐马尔科夫模型，一共生成10个隐马尔科夫模型，其中每个隐马尔科夫模型包含三个参数；
(6.3)、将每个随机生成的隐马尔科夫模型分别导入公式（13），得到每个隐马尔科夫模型在公式（13）中的值，找到全局最好位置，并记录全局最好位置的适应值；利用粒子移动公式进行粒子的更新，公式如下：
Figure DEST_PATH_IMAGE118
(14)
Figure DEST_PATH_IMAGE120
(15)
Figure DEST_PATH_IMAGE122
（16）
式中
Figure DEST_PATH_IMAGE124
—第i个粒子在第t次迭代的最好位置矢量；
Figure DEST_PATH_IMAGE126
—第t次迭代时的全局最好位置矢量
Figure DEST_PATH_IMAGE128
—第t次迭代时全体粒子当前最好位置的中心位置；
Figure DEST_PATH_IMAGE130
与
Figure 298352DEST_PATH_IMAGE124
之间的随机位置；
Figure DEST_PATH_IMAGE132
和
Figure DEST_PATH_IMAGE134
为[0,1]上均匀分布的随机数；
Figure DEST_PATH_IMAGE136
是收缩—扩张系数，它是算法中群体规模和迭代次数以外的唯一一个可控参数，对于不同取值的
Figure 341132DEST_PATH_IMAGE136
，会对粒子的收敛性产生影响；通常采取线性减少的方式从m减少到n，其公式如下。
Figure DEST_PATH_IMAGE138
（17）
一般取m=1，n=0.5，
Figure DEST_PATH_IMAGE140
为最大迭代次数；
(6.4)、当每次迭代得到的适应值保持一致时，迭代结束，当前适应值最大的粒子为隐马尔科夫模型训练的最优化参数。
所述的基于HMM和QPSO优化算法的轴承故障程度辨识方法，所述诊断轴承，包括(7.1)、用上述步骤(2)至(5)对测试信号进行处理，得到测试信号的观测数列；
(7.2)、用前向算法输出的最大似然估计确定每组测试信号所对应的故障程度；其算法流程如下。
Figure DEST_PATH_IMAGE142
（18）
Figure DEST_PATH_IMAGE144
（19）
Figure DEST_PATH_IMAGE146
（20）
(7.3)、根据每组测试信号在各个隐马尔可夫模型中得到的最大似然值，其最大似然值越大，代表其和当前隐马尔可夫模型对应的故障程度越相近。
本发明的优点与效果是
本发明用参数作为隐马尔科夫模型的学习参数，从而对轴承故障程度进行分类，用以故障程度评估，具有重要的工程意义。
附图说明
图1为本发明提出的故障诊断模型流程图；
图2为本发明重度故障信号变分模态分解分解后的IMF示意图；
图3为本发明奇异值分解之后部分奇异值表；
图4为本发明不同种类的轴承故障信号图；
图5为本发明轴承数据采集设备图；
图6 为本发明 QPSO算法迭代图；
图7为本发明正常隐马尔科夫模型下，不同损坏程度的轴承最大似然；
图8为本发明轻度故障隐马尔科夫模型下，不同损坏程度的轴承最大似然；
图9为本发明中度故障隐马尔科夫模型下，不同损坏程度的轴承最大似然；
图10为本发明重度故障隐马尔科夫模型下，不同损坏程度的轴承最大似然。
具体实施方式
下面结合附图所示实施例对本发明进行详细说明。
总体的，如图1所示，本发明实施例公开了一种基于HMM和QPSO优化算法的轴承故障程度辨识方法，包括如下步骤：
1)、采集滚动轴承振动信号使用振动数据采集仪以采样频率为12000Hz采集待检测滚动轴承在运行状态下的振动信号，记为X[m]，m＝1 ,2 ,…,M，M为总采样点数，并标记轴承状态，总状态数为y，将所有的所得的数据集合为
Figure DEST_PATH_IMAGE148
，标签数据集为
Figure DEST_PATH_IMAGE150
。
(2)、数据的分帧操作，建立训练集
取前120k个采样点，平均分成12组数据，每组数据10k个采样点，再把每组数据平均分成5帧，每帧2k个采样点。
(3)、振动信号的分频处理
利用变分模态分解VMD的可变尺度分解对训练集中的信号进行处理，变分模态分解的过程是将一个复杂的滚动轴承时域振动信号分解为K个IMF分量，K为任意正整数，变分模态分解的基本过程如下：
(3.1)、希尔伯特变换，对每一个IMF分量进行希尔伯特变换，得到单边频谱如下：
Figure RE-DEST_PATH_IMAGE081
（1）
式中
Figure RE-RE-DEST_PATH_IMAGE082
是第k个IMF分量，k为任意正整数，
Figure RE-365553DEST_PATH_IMAGE082
可看做谐波信号：
Figure RE-DEST_PATH_IMAGE083
，其中
Figure RE-RE-DEST_PATH_IMAGE084
是
Figure RE-879884DEST_PATH_IMAGE082
的瞬时幅值，且
Figure RE-DEST_PATH_IMAGE085
为
Figure RE-30374DEST_PATH_IMAGE082
的相位，且
Figure RE-RE-DEST_PATH_IMAGE086
为冲激函数；
Figure RE-DEST_PATH_IMAGE087
；t是时间；
Figure RE-RE-DEST_PATH_IMAGE088
为卷积运算。
(3.2)、频率混合，利用上述得出的单边频谱信号混合一个估计的中心频率
Figure RE-DEST_PATH_IMAGE089
，待混合之后，每一个
Figure RE-RE-DEST_PATH_IMAGE090
的频谱将移动到相应的基频带上，即：
Figure RE-DEST_PATH_IMAGE091
（2）
(3.3)、带宽估计，计算公式（2）所得到的解调信号的梯度二范数平方：
Figure RE-RE-DEST_PATH_IMAGE092
（3）
其中
Figure RE-DEST_PATH_IMAGE093
为函数对时间t的偏导数；
(3.4)、建立最优化模型，通过以上计算，可以得到的变分约束模型如下：
Figure RE-RE-DEST_PATH_IMAGE094
（4）
式中，K为IMF分量个数；
Figure RE-210466DEST_PATH_IMAGE090
是模态的集合，
Figure RE-DEST_PATH_IMAGE095
是
Figure RE-420998DEST_PATH_IMAGE090
的中心频率；
Figure RE-RE-DEST_PATH_IMAGE096
是K个IMF
Figure RE-25286DEST_PATH_IMAGE090
的加和；
(3.5)、转化为非约束问题，引入二次惩罚因子
Figure RE-DEST_PATH_IMAGE097
和拉格朗日乘子
Figure RE-RE-DEST_PATH_IMAGE098
构成表达式：
Figure RE-DEST_PATH_IMAGE099
（5）
Figure RE-RE-DEST_PATH_IMAGE100
，
Figure RE-DEST_PATH_IMAGE101
和
Figure RE-RE-DEST_PATH_IMAGE102
初始值都为0，通过交替方向乘子法（ADMM）迭代得到（4）的最优解，也可以得到
Figure RE-480932DEST_PATH_IMAGE100
所对应的中心频率
Figure RE-191678DEST_PATH_IMAGE101
如下所示：
Figure RE-DEST_PATH_IMAGE103
（6）
Figure RE-RE-DEST_PATH_IMAGE104
（7）
式中：
Figure RE-DEST_PATH_IMAGE105
为当前余项
Figure RE-RE-DEST_PATH_IMAGE106
的维纳滤波，
Figure RE-DEST_PATH_IMAGE107
为模态功率谱的重心；
Figure RE-RE-DEST_PATH_IMAGE108
为
Figure RE-DEST_PATH_IMAGE109
的傅里叶变换，n为任意正整数，
Figure RE-RE-DEST_PATH_IMAGE110
为频率，i为任意正整数；
(3.6)、变分模态分解（VMD）对滚动轴承轴承振动信号分解过程如下：设置n=n+1，k=k+1，依照公式（6）（7）更新
Figure RE-957465DEST_PATH_IMAGE100
和
Figure RE-908104DEST_PATH_IMAGE101
，直到k=K是结束迭代；依照公式（8）更新
Figure RE-DEST_PATH_IMAGE111
公式如下：
Figure RE-RE-DEST_PATH_IMAGE112
（8）
式中
Figure RE-DEST_PATH_IMAGE113
为噪声容限参数，通常令
Figure RE-RE-DEST_PATH_IMAGE114
，重复上述操作直到满足公式（9），就可以得到由变分模态分解（VMD）处理后的最优化IMF；如图2所示。公式如下：
Figure RE-DEST_PATH_IMAGE115
（9）
式中
Figure RE-RE-DEST_PATH_IMAGE116
为误差，可以取任意正整数。
(4)、对本征模态分量进行奇异值分解（SVD）
利用奇异值分解（SVD），提取由变分模态分解（VMD）分解后的最优IMF的奇异值矩阵如图3所示。奇异值分解公式如下：
Figure RE-DEST_PATH_IMAGE117
（10）
其中
Figure RE-RE-DEST_PATH_IMAGE118
为左奇异向量矩阵，
Figure RE-DEST_PATH_IMAGE119
为右奇异向量矩阵，
Figure RE-RE-DEST_PATH_IMAGE120
为奇异值向量矩阵，它是一个对角矩阵，其主对角线上的值为矩阵IMF的奇异值，其他位置上的元素值都为0。用此方法，分别求出每组数据的奇异值向量。
(5)、对奇异值向量进行聚类
使用K-means聚类算法，将每组数据进行分类，其算法流程如下：(5.1)、从数据集中随机选取k个样本作为初始的k个质心向量：
Figure DEST_PATH_IMAGE160
；(5.2)、对于
Figure DEST_PATH_IMAGE162
(5.2.1)、将簇划分
Figure DEST_PATH_IMAGE164
初始化为
Figure DEST_PATH_IMAGE166
，
Figure DEST_PATH_IMAGE168
；
(5.2.2)、对于
Figure DEST_PATH_IMAGE170
，计算样本
Figure DEST_PATH_IMAGE172
和各个质心向量
Figure DEST_PATH_IMAGE174
的距离，公式如下：
Figure DEST_PATH_IMAGE176
（11）
将
Figure DEST_PATH_IMAGE178
标记最小的
Figure DEST_PATH_IMAGE180
所对应的类别
Figure DEST_PATH_IMAGE182
，此时更新
Figure DEST_PATH_IMAGE184
；
(5.2.3)、对于
Figure DEST_PATH_IMAGE186
，对
Figure DEST_PATH_IMAGE188
中所有的样本点重新计算新的质心，公式如下：
Figure DEST_PATH_IMAGE190
（12）
(5.2.4)、如果所有的k个质心向量都没有发生变化则转到步骤（5.2.3）
(5.3)输出簇划分；由此得到了一组整数，作为隐马尔可夫模型的观测数据，用于后期的故障辨识。
(6)、利用量子粒子群算法对隐马尔科夫模型的参数
Figure DEST_PATH_IMAGE192
进行估计。每一个粒子都是一个HMM参数，令其适应度函数为：
Figure DEST_PATH_IMAGE194
(13)
利用量子粒子群算法求隐马尔科夫模型的最优化参数步骤如下：(6.1)、初始化，种群数量10，最大迭代次数30，算法运行10轮； (6.2)、随机生成隐马尔科夫模型，一共生成10个隐马尔科夫模型，其中每个隐马尔科夫模型包含
Figure 309045DEST_PATH_IMAGE192
三个参数。
(6.3)、将每个随机生成的隐马尔科夫模型分别导入公式（13），得到每个隐马尔科夫模型在公式（13）中的值，找到全局最好位置，并记录全局最好位置的适应值。利用粒子移动公式进行粒子的更新，公式如下：
Figure DEST_PATH_IMAGE196
(14)
Figure DEST_PATH_IMAGE198
(15)
Figure DEST_PATH_IMAGE200
(16)
式中
Figure DEST_PATH_IMAGE202
—第i个粒子在第t次迭代的最好位置矢量；
Figure 269304DEST_PATH_IMAGE126
—第t次迭代时的全局最好位置矢量
Figure 332069DEST_PATH_IMAGE128
—第t次迭代时全体粒子当前最好位置的中心位置；
Figure 317342DEST_PATH_IMAGE130
与
Figure 689418DEST_PATH_IMAGE124
之间的随机位置；
Figure 71727DEST_PATH_IMAGE132
和
Figure 125133DEST_PATH_IMAGE134
为[0,1]上均匀分布的随机数；
Figure 964913DEST_PATH_IMAGE136
是收缩—扩张系数，它是算法中群体规模和迭代次数以外的唯一一个可控参数，对于不同取值的
Figure 507890DEST_PATH_IMAGE136
，会对粒子的收敛性产生影响；通常采取线性减少的方式从m减少到n，其公式如下。
Figure 800331DEST_PATH_IMAGE138
（17）
一般取m=1，n=0.5，
Figure 847309DEST_PATH_IMAGE140
为最大迭代次数；
(6.4)、当每次迭代得到的适应值保持一致时，如图6所示，迭代结束，当前适应值最大的粒子为隐马尔科夫模型训练的最优化参数。
(7)、诊断轴承
7.1)、用上述步骤(2)至(5)对测试信号进行处理，得到测试信号的观测数列。
(7.2)、用前向算法输出的最大似然估计确定每组测试信号所对应的故障程度。其算法流程如下。
Figure 541596DEST_PATH_IMAGE142
（18）
Figure 255474DEST_PATH_IMAGE144
（19）
Figure 113839DEST_PATH_IMAGE146
（20）
(7.3)、根据每组测试信号在各个隐马尔可夫模型中得到的最大似然值，其最大似然值越大，代表其和当前隐马尔可夫模型对应的故障程度越相近，如图7、8、9、10所示。
(8)、诊断结束。
在实施例中，依托试验室中的轴承振动试验台，试验台如图5 所示。
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (4)

1.基于HMM和QPSO优化算法的轴承故障程度辨识方法，其特征在于，所述方法包括如下步骤：
(1)、采集滚动轴承振动信号
使用振动数据采集仪以采样频率为12000Hz采集待检测滚动轴承在运行状态下的振动信号，记为X[m]，m＝1 ,2 ,…,M，M为总采样点数，并标记轴承状态，总状态数为y，将所有的所得的数据集合为
Figure 483070DEST_PATH_IMAGE001
，标签数据集为
Figure 237400DEST_PATH_IMAGE002
；
(2)、数据的分帧操作，建立训练集
取前120k个采样点，平均分成12组数据，每组数据10k个采样点，再把每组数据平均分成5帧，每帧2k个采样点；
(3)、振动信号的分频处理
利用变分模态分解VMD的可变尺度分解对训练集中的信号进行处理，变分模态分解的过程是将一个复杂的滚动轴承时域振动信号分解为K个IMF分量，K为任意正整数；
(4)、对本征模态分量进行奇异值分解（SVD）；利用奇异值分解（SVD），提取由变分模态分解（VMD）分解后的最优IMF的奇异值矩阵；
(5)、对奇异值向量进行聚类；
(6)、利用量子粒子群算法对隐马尔科夫模型的参数
Figure DEST_PATH_IMAGE003
进行估计；每一个粒子都是一个HMM参数，令其适应度函数为：
Figure 36729DEST_PATH_IMAGE004
（13）
(7)、诊断轴承；
(8)、诊断结束；
所述利用量子粒子群算法求隐马尔科夫模型的最优化参数步骤如下：
(6.1)、初始化，种群数量10，最大迭代次数30，算法运行10轮；
(6.2)、随机生成隐马尔科夫模型，一共生成10个隐马尔科夫模型，其中每个隐马尔科夫模型包含
Figure DEST_PATH_IMAGE005
三个参数；
(6.3)、将每个随机生成的隐马尔科夫模型分别导入公式（13），得到每个隐马尔科夫模型在公式（13）中的值，找到全局最好位置，并记录全局最好位置的适应值；利用粒子移动公式进行粒子的更新，公式如下：
Figure 622431DEST_PATH_IMAGE006
（14）
Figure DEST_PATH_IMAGE007
（15）
Figure 616057DEST_PATH_IMAGE008
（16）
式中
Figure 857682DEST_PATH_IMAGE009
—第i个粒子在第t次迭代的最好位置矢量；
Figure DEST_PATH_IMAGE010
—第t次迭代时的全局最好位置矢量
Figure 257440DEST_PATH_IMAGE011
—第t次迭代时全体粒子当前最好位置的中心位置；
Figure DEST_PATH_IMAGE012
与
Figure 963227DEST_PATH_IMAGE013
之间的随机位置；
Figure DEST_PATH_IMAGE014
和
Figure 865105DEST_PATH_IMAGE015
为[0,1]上均匀分布的随机数；
Figure DEST_PATH_IMAGE016
是收缩—扩张系数，它是算法中群体规模和迭代次数以外的唯一一个可控参数，对于不同取值的
Figure 656344DEST_PATH_IMAGE016
，会对粒子的收敛性产生影响；通常采取线性减少的方式从m减少到n，其公式如下：
Figure 735158DEST_PATH_IMAGE017
（17）
一般取m=1，n=0.5，
Figure DEST_PATH_IMAGE018
为最大迭代次数；
(6.4)、当每次迭代得到的适应值保持一致时，迭代结束，当前适应值最大的粒子为隐马尔科夫模型训练的最优化参数；
所述诊断轴承，包括(7.1)、用上述步骤(2)至(5)对测试信号进行处理，得到测试信号的观测数列；
(7.2)、用前向算法输出的最大似然估计确定每组测试信号所对应的故障程度；其算法流程如下：
Figure 295453DEST_PATH_IMAGE019
（18）
Figure DEST_PATH_IMAGE020
（19）
Figure 129416DEST_PATH_IMAGE021
(20)
(7.3)、根据每组测试信号在各个隐马尔可夫模型中得到的最大似然值，其最大似然值越大，代表其和当前隐马尔可夫模型对应的故障程度越相近。
2.根据权利要求1所述的基于HMM和QPSO优化算法的轴承故障程度辨识方法，其特征在于，所述变分模态分解的基本过程如下：
(3.1)、希尔伯特变换，对每一个IMF分量进行希尔伯特变换，得到单边频谱如下：
Figure DEST_PATH_IMAGE022
(1)
式中
Figure 643837DEST_PATH_IMAGE023
是第k个IMF分量，k为任意正整数，
Figure 588659DEST_PATH_IMAGE023
看做谐波信号：
Figure DEST_PATH_IMAGE024
，其中
Figure 3460DEST_PATH_IMAGE025
是
Figure 946008DEST_PATH_IMAGE023
的瞬时幅值，且
Figure DEST_PATH_IMAGE026
为
Figure 446259DEST_PATH_IMAGE023
的相位，且
Figure 427728DEST_PATH_IMAGE027
；
Figure DEST_PATH_IMAGE028
为冲激函数；
Figure 962615DEST_PATH_IMAGE029
；t是时间；
Figure DEST_PATH_IMAGE030
为卷积运算；
(3.2)、频率混合，利用上述得出的单边频谱信号混合一个估计的中心频率
Figure 138381DEST_PATH_IMAGE031
，待混合之后，每一个
Figure 63612DEST_PATH_IMAGE023
的频谱将移动到相应的基频带上，即：
Figure DEST_PATH_IMAGE032
(2)
(3.3)、带宽估计，计算公式（2）所得到的解调信号的梯度二范数平方：
Figure 350237DEST_PATH_IMAGE033
（3）
其中
Figure DEST_PATH_IMAGE034
为函数对时间t的偏导数；
(3.4)、建立最优化模型，通过以上计算，可以得到的变分约束模型如下：
Figure 975516DEST_PATH_IMAGE035
(4)
式中，K为IMF分量个数；
Figure DEST_PATH_IMAGE036
是模态的集合，
Figure 322183DEST_PATH_IMAGE037
是
Figure 734710DEST_PATH_IMAGE036
的中心频率；
Figure DEST_PATH_IMAGE038
是K个IMF
Figure 825026DEST_PATH_IMAGE036
的加和；
(3.5)、转化为非约束问题，引入二次惩罚因子
Figure 68925DEST_PATH_IMAGE039
和拉格朗日乘子
Figure DEST_PATH_IMAGE040
构成表达式：
Figure 813591DEST_PATH_IMAGE041
（5）
Figure DEST_PATH_IMAGE042
，
Figure 775731DEST_PATH_IMAGE043
和
Figure DEST_PATH_IMAGE044
初始值都为0，通过交替方向乘子法（ADMM）迭代得到（4）的最优解，也可以得到
Figure 404158DEST_PATH_IMAGE042
所对应的中心频率
Figure 705826DEST_PATH_IMAGE043
如下所示：
Figure 128718DEST_PATH_IMAGE045
（6）
Figure DEST_PATH_IMAGE046
（7）
式中：
Figure 79618DEST_PATH_IMAGE047
为当前余项
Figure DEST_PATH_IMAGE048
的维纳滤波，
Figure 511737DEST_PATH_IMAGE049
为模态功率谱的重心；
Figure DEST_PATH_IMAGE050
为
Figure 464649DEST_PATH_IMAGE051
的傅里叶变换，n为任意正整数，
Figure DEST_PATH_IMAGE052
为频率，i为任意正整数；
(3.6)、变分模态分解（VMD）对滚动轴承轴承振动信号分解过程如下：设置n=n+1，k=k+1，
依照公式（6）（7）更新
Figure 324021DEST_PATH_IMAGE053
和
Figure DEST_PATH_IMAGE054
，直到k=K是结束迭代；依照公式（8）更新
Figure 493709DEST_PATH_IMAGE055
公式如下：
Figure DEST_PATH_IMAGE056
（8）
式中
Figure 729518DEST_PATH_IMAGE057
为噪声容限参数，通常令
Figure DEST_PATH_IMAGE058
，重复上述操作直到满足公式（9），就可以得到由变分模态分解（VMD）处理后的最优化IMF；公式如下：
Figure 802516DEST_PATH_IMAGE059
（9）
式中
Figure DEST_PATH_IMAGE060
为误差，可以取任意正整数。
3.根据权利要求1所述的基于HMM和QPSO优化算法的轴承故障程度辨识方法，其特征在于，所述奇异值分解公式如下：
Figure 832789DEST_PATH_IMAGE061
（10）
其中
Figure DEST_PATH_IMAGE062
为左奇异向量矩阵，
Figure 492703DEST_PATH_IMAGE063
为右奇异向量矩阵，
Figure DEST_PATH_IMAGE064
为奇异值向量矩阵，它是一个对角矩阵，其主对角线上的值为矩阵IMF的奇异值，其他位置上的元素值都为0；用此方法，分别求出每组数据的奇异值向量。
4.根据权利要求1所述的基于HMM和QPSO优化算法的轴承故障程度辨识方法，其特征在于，所述对奇异值向量进行聚类，使用K-means聚类算法，将每组数据进行分类，其算法流程如下：
（1）从数据集中随机选取k个样本作为初始的k个质心向量：
Figure 266624DEST_PATH_IMAGE065
；
（2）对于
Figure DEST_PATH_IMAGE066
（a）将簇划分
Figure 131812DEST_PATH_IMAGE067
初始化为
Figure DEST_PATH_IMAGE068
，
Figure 332986DEST_PATH_IMAGE069
Figure DEST_PATH_IMAGE070
；
（b）对于
Figure 978731DEST_PATH_IMAGE071
，计算样本
Figure DEST_PATH_IMAGE072
和各个质心向量
Figure 795158DEST_PATH_IMAGE073
的距离，公式如下：
Figure DEST_PATH_IMAGE074
（11）
将
Figure 842748DEST_PATH_IMAGE075
标记最小的
Figure DEST_PATH_IMAGE076
所对应的类别
Figure 214824DEST_PATH_IMAGE077
，此时更新
Figure DEST_PATH_IMAGE078
；
（c）对于
Figure 347865DEST_PATH_IMAGE079
，对
Figure DEST_PATH_IMAGE080
中所有的样本点重新计算新的质心，公式如下：
Figure 837490DEST_PATH_IMAGE081
（12）
（d）如果所有的k个质心向量都没有发生变化则转到步骤（c）
（3）输出簇划分；由此得到了一组整数，作为隐马尔可夫模型的观测数据，用于后期的故障辨识。
CN201911342148.2A 2019-12-24 2019-12-24 基于hmm和qpso优化算法的轴承故障程度辨识方法 Active CN111189638B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911342148.2A CN111189638B (zh) 2019-12-24 2019-12-24 基于hmm和qpso优化算法的轴承故障程度辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911342148.2A CN111189638B (zh) 2019-12-24 2019-12-24 基于hmm和qpso优化算法的轴承故障程度辨识方法

Publications (2)

Publication Number Publication Date
CN111189638A CN111189638A (zh) 2020-05-22
CN111189638B true CN111189638B (zh) 2021-08-06

Family

ID=70707434

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911342148.2A Active CN111189638B (zh) 2019-12-24 2019-12-24 基于hmm和qpso优化算法的轴承故障程度辨识方法

Country Status (1)

Country Link
CN (1) CN111189638B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111753927A (zh) * 2020-07-22 2020-10-09 华东交通大学 一种机车牵引座裂纹故障状态识别方法
CN112014047B (zh) * 2020-08-27 2022-05-03 华侨大学 一种有载分接开关机械故障诊断方法
CN113340625B (zh) * 2021-04-21 2023-01-06 北京交通大学 转向架故障诊断方法
CN114252261B (zh) * 2021-12-21 2024-03-08 沈阳顺义科技股份有限公司 一种综合传动装置转向系统故障诊断方法及系统

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104198924A (zh) * 2014-09-11 2014-12-10 合肥工业大学 一种新颖的模拟电路早期故障诊断方法
CN104361414A (zh) * 2014-11-24 2015-02-18 武汉大学 一种基于相关向量机的输电线路覆冰预测方法
CN105258839A (zh) * 2015-10-30 2016-01-20 南京信息工程大学 基于量子粒子群小波神经网络的阵列式气压测量补偿装置及其方法
CN108320516A (zh) * 2018-04-08 2018-07-24 华中师范大学 基于尖点突变和量子粒子群优化的道路通行能力评价方法
CN108898602A (zh) * 2018-06-27 2018-11-27 南京邮电大学 一种基于改进qpso的fcm医学图像分割方法
CN109029975A (zh) * 2018-06-26 2018-12-18 红河学院 一种风电齿轮箱的故障诊断方法
CN109580222A (zh) * 2018-12-04 2019-04-05 河北科技大学 基于变分模态分解-传递熵的轴承退化状态识别预测方法
CN109992844A (zh) * 2019-03-13 2019-07-09 上海电力学院 一种基于adqpso-svr模型的锅炉飞灰含碳量预测方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102054179A (zh) * 2010-12-14 2011-05-11 广州大学 一种旋转机械在线状态监测与故障诊断装置及方法
KR101877141B1 (ko) * 2016-12-08 2018-07-10 한국항공우주연구원 안테나 패턴 합성 장치 및 방법
CN107378641B (zh) * 2017-08-23 2019-02-01 东北电力大学 一种基于图像特征和lltsa算法的刀具磨损状态监测方法
CN108171263B (zh) * 2017-12-26 2019-08-30 合肥工业大学 基于改进变分模态分解和极限学习机的滚动轴承故障诊断方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104198924A (zh) * 2014-09-11 2014-12-10 合肥工业大学 一种新颖的模拟电路早期故障诊断方法
CN104361414A (zh) * 2014-11-24 2015-02-18 武汉大学 一种基于相关向量机的输电线路覆冰预测方法
CN105258839A (zh) * 2015-10-30 2016-01-20 南京信息工程大学 基于量子粒子群小波神经网络的阵列式气压测量补偿装置及其方法
CN108320516A (zh) * 2018-04-08 2018-07-24 华中师范大学 基于尖点突变和量子粒子群优化的道路通行能力评价方法
CN109029975A (zh) * 2018-06-26 2018-12-18 红河学院 一种风电齿轮箱的故障诊断方法
CN108898602A (zh) * 2018-06-27 2018-11-27 南京邮电大学 一种基于改进qpso的fcm医学图像分割方法
CN109580222A (zh) * 2018-12-04 2019-04-05 河北科技大学 基于变分模态分解-传递熵的轴承退化状态识别预测方法
CN109992844A (zh) * 2019-03-13 2019-07-09 上海电力学院 一种基于adqpso-svr模型的锅炉飞灰含碳量预测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Multiple sequence alignment using the Hidden Markov Model trained by an improved quantum-behaved particle swarm optimization;Jun Sun;《Information sciences》;20121231;全文 *
一种基于QPSO_RVM的模拟电路故障预测方法;张朝龙;《仪器仪表学报》;20140831;全文 *
基于量子粒子群优化算法的机器人运动学标定方法;房立金;《机械工程学报》;20160430;全文 *

Also Published As

Publication number Publication date
CN111189638A (zh) 2020-05-22

Similar Documents

Publication Publication Date Title
CN111189638B (zh) 基于hmm和qpso优化算法的轴承故障程度辨识方法
US20200271720A1 (en) Method for diagnosing analog circuit fault based on vector-valued regularized kernel function approximation
CN103728551B (zh) 一种基于级联集成分类器的模拟电路故障诊断方法
CN102629374B (zh) 基于子空间投影和邻域嵌入的图像超分辨率重建方法
CN108171119B (zh) 基于残差网络的sar图像变化检测方法
CN106384092A (zh) 面向监控场景的在线低秩异常视频事件检测方法
CN111238843B (zh) 一种基于快速谱峭度分析的风机健康评价方法
CN111476339B (zh) 滚动轴承故障特征提取方法、智能诊断方法及系统
CN110020714B (zh) 模型训练及数据分析方法、装置、设备以及存储介质
CN114282571B (zh) 一种轴承多维健康指标构建方法、系统、设备以及介质
CN112287796B (zh) 基于VMD-Teager能量算子的辐射源识别方法
CN110738232A (zh) 一种基于数据挖掘技术的电网电压越限成因诊断方法
CN111310719B (zh) 一种未知辐射源个体识别及检测的方法
CN114429152A (zh) 基于动态指数对抗性自适应的滚动轴承故障诊断方法
CN114169377A (zh) 基于g-mscnn的有噪环境中滚动轴承故障诊断方法
CN113920255B (zh) 基于点云数据的高效测绘系统
CN111665050A (zh) 一种基于聚类k-svd算法的滚动轴承故障诊断方法
CN110929761A (zh) 智能系统安全体系态势感知构架中采集样本的平衡方法
Deng et al. Impulse feature extraction method for machinery fault detection using fusion sparse coding and online dictionary learning
CN110718301A (zh) 基于动态脑功能网络的阿尔茨海默病辅助诊断装置及方法
CN113758709A (zh) 结合边缘计算和深度学习的滚动轴承故障诊断方法及系统
Guo et al. An equipment multiple failure causes intelligent identification method based on integrated strategy for subway sliding plug door system under variable working condition
CN112464712B (zh) 一种基于盲抽取算法的旋转机械故障诊断方法
CN104408072B (zh) 一种基于复杂网络理论的适用于分类的时间序列特征提取方法
Li et al. Fault diagnosis of ZD6 turnout system based on wavelet transform and GAPSO-FCM

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant