CN107144428A - 一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法 - Google Patents

一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法 Download PDF

Info

Publication number
CN107144428A
CN107144428A CN201710160858.8A CN201710160858A CN107144428A CN 107144428 A CN107144428 A CN 107144428A CN 201710160858 A CN201710160858 A CN 201710160858A CN 107144428 A CN107144428 A CN 107144428A
Authority
CN
China
Prior art keywords
mrow
msub
mtr
mtd
value
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.)
Granted
Application number
CN201710160858.8A
Other languages
English (en)
Other versions
CN107144428B (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.)
Beijing Jiaotong University
Original Assignee
Beijing Jiaotong University
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 Beijing Jiaotong University filed Critical Beijing Jiaotong University
Priority to CN201710160858.8A priority Critical patent/CN107144428B/zh
Publication of CN107144428A publication Critical patent/CN107144428A/zh
Application granted granted Critical
Publication of CN107144428B publication Critical patent/CN107144428B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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

Abstract

本发明公开一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，包括：S100、进行多角度特征值提取与融合；S200、进行基于循环平稳理论的故障诊断；S300、通过进行多变量模型训练实现基于多变量模型的故障分离；S400、基于半监督算法，协同使用BP神经网络和支持向量回归两种算法进行基于故障诊断的剩余寿命预测。本发明实现了轴承发生早期故障的故障检测、故障分离及剩余寿命预测，为车辆轴承维修决策的制定提供依据。

Description

一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法
技术领域
本发明涉及轨道交通车辆轴承的故障诊断及剩余寿命预测领域。更具体地，涉及一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法。
背景技术
轨道交通的高速发展，对车辆的安全性和可靠性提出了更高的要求。早晚高峰时段很多运营线路均在超负荷运行，这就要求列车各个关键系统必须安全可靠，不能出现故障或是失效状态下运行的现象。轴承是轨道交通中不可缺少的元件之一，同时也是最易损坏的元件，其有效性直接导致列车的运行安全。复杂时变的运行环境：载荷，轨道的平顺程度、温度等，对于轴承的健康状态有很大的影响。当轴承中某一部位发生故障时，会产生连锁反应，轻则导致车辆的相关系统受损，重则会导致车辆的停运甚至是人员的伤亡。
在现有的剩余寿命预测研究中，较少的考虑到故障类型对于剩余寿命预测的影响。当故障发生后轴承进入快速衰退状态，而不同故障状态下轴承的衰退特性大不相同，剩余寿命预测的结果也大不相同。
因此，需要提供一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法。
发明内容
本发明的目的在于提供一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，以实现轴承发生早期故障的故障检测、故障分离及剩余寿命预测，为车辆轴承维修决策的制定提供依据。
为达到上述目的，本发明采用下述技术方案：
一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，包括如下步骤：
S100、使用局部均值分解法将轴承振动信号分解为若干PF单分量，对分离出的单分量进行时域、频域、能量及统计量的多角度特征值提取，将提取到的多角度特征值使用主成分分析法降维与融合；
S200、进行基于循环平稳理论的故障诊断；
S300、通过进行多变量模型训练实现基于多变量模型的故障分离；
S400、基于半监督算法，协同使用BP神经网络和支持向量回归两种算法进行基于故障诊断的剩余寿命预测。
优选地，所述使用局部均值分解法将轴承振动信号分解为若干PF单分量进一步包括：
S111、寻找轴承振动信号x(t)所有局部极值点ni，将上下极值点分别用三次样条曲线进行连接，得到信号的上、下包络线Emax、Emin，计算局部均值函数m11(t)和包络函数c11(t)：
S112、从x(t)中分离出局部均值函数m11(t)，得到：
h11(t)＝x(t)-m11(t)；
S113、对h11(t)解调，将h11(t)除以包络函数c11(t)得到：
S114、按照S111中的方法求出s11(t)所对应的包络函数c12(t)，若包络函数c12(t)＝1，则判断s11(t)为纯调频函数；若包络函数c12(t)≠1，则迭代执行上述步骤S111-S113，直至s1n(t)的包络估计函数c1(n+1)(t)＝1，则有：
S115、将进行迭代处理过程中产生的所有包络函数相乘得包络信号：
S116、令包络信号c1(t)与纯调频函数s1n(t)相乘，得到x(t)的第一个PF分量PF1：
PF1(t)＝c1(t)·s1n(t)；
S117、从x(t)中分离出第一个分量PF1后得到一个新的待分解信号r1(t)，将新的待分解信号替代轴承振动信号迭代执行步骤S111-S116，直至新的待分解信号rk(t)为一单调函数为止，k代表迭代次数，即
至此，x(t)被分解为k个PF分量及一个单调函数rk(t)之和：
优选地，所述多角度特征值包括时域特征值、频域特征值、统计特征参数特征值和能量特征值；所述时域特征值包括均方值、方差、峰值、峭度、偏度、脉冲因子、峰值因子和裕度因子；所述频域特征值包括均方频率、频谱重心、频率方差和频率标准差；所述统计特征参数特征值包括威布尔形状参数、威布尔尺度参数、伽马形状参数和伽马尺度参数；所述能量特征值包括香农熵和能量矩。
优选地，所述将提取到的多角度特征值使用主成分分析法降维与融合进一步包括：
S131、设多角度特征值构成的数据的样本数为N，每个样本包含P个特征向量，样本矩阵为XNxP，对数据进行标准化：
其中，i＝1,2,…,n；；j＝1,2,…,p；；
S132、计算特征向量的协方差矩阵R，并计算其特征值和特征向量：
R＝YTY/(N-1)
其中，Y是标准化后的数据，计算出R的特征值分别为λ1≥λ2≥…≥λp，以及其对应的特征向量为αi＝(αi1,αi2,…,αip)T,i＝1,2,…p；
S133、设正交空间中前k个主元分量为y1,y2…yk，计算累计方差贡献率：
至此，完成了将P维数据降至K维数据，实现了多角度特征值降维与融合。
优选地，步骤S200进一步包括：
S211、对于为非平稳随机信号的轴承振动信号x(t)，x(t)的时变自相关函数的表示为：
其中，τ为时延因子，E{·}表示统计平均，*表示复共轭；
S212、设Rx(t,τ)以T为周期，用样本平均代替统计平均，将时变自相关函数表示为：
将上式用傅里叶级数展开表示为：
其中α＝m/T(m∈Z)为循环频率，其傅里叶系数为
S213、设T0＝(2N+1)T，通过整理可以得到循环自相关函数为：
R(τ,α)＝<x(t+τ/2)x*(t-τ/2)e-j2παt>t；
S214、得到循环谱相关密度函数并根据循环谱相关密度函数进行基于循环平稳理论的故障诊断，循环谱相关密度函数为：
其中，f为谱频率。
优选地，步骤S300进一步包括：
进行多变量模型训练：
S311、针对具有H种故障类型的数据，样本总数为N，每种故障下样本数为[n1,n2,…nH]；
S312、对于不同故障状态的所有样本进行振动数据特征提取，使用拉普拉斯分值法对所有特征值进行选取，每种故障都选出p个最优特征量，分别为X＝[X1,X2,…,Xp]
S313、设模型类型m＝1，0＜m＜4，模型阶数r＝1，0＜r＜p-1，设h＝1；
S314、对第h(1≤h≤H)类训练样本执行：选择特征量Xi为被预测变量，选择r个特征量Xj(j≠i)为预测变量，则共有种选择方案，即对于被预测变量来说可建立中预测模型；
S315、第h类故障类型共有nh个训练样本，对于每种组合方式，被预测变量Xi都可以建立nh个方程，使用最小二乘法对模型bo,bj,bjj,bjk进行参数估计，得到每种组合下的被预测变量的数学方程
S316、计算每种组合下所有样本对Xi的估计误差和其中l表示第l个训练样本，共有个SSE值，选择最小的SSE对应的预测模型为该种故障，则模型类型为m且阶数为r的条件下，变量Xi的预测模型已确定；
S317、令h＝h+1，迭代执行步骤S313-S316，直到h＝H；
S318、令r＝r+1，迭代执行S314-S317直到r＝P-1；
S319、令m＝m+1，迭代执行S314-S318直到m＝4，得到了各种模型类型及各种阶次条件下的预测模型即每个m和r的条件下都会得到H*P个组成矩阵，然后将训练变量回代各个方程，选择最小误差值下的m和r对应的矩阵为本次的训练输出模型；
进行故障类型预测：
S321、对预测样本计算同样的p个最优特征量，分别为X＝[X1,X2,…,Xp]；
S322、将特征值带入训练好的数学模型，得到预测的估计变量值其中h＝[1,2,…,H]，i＝[1,2,…,p]；
S323、计算每种类型下，所有估计变量的误差值之和SSEh最小误差值的状态类型为该待预测数据的故障类型。
优选地，步骤S400进一步包括：
S411、将故障数据等间隔取值，3/4作为训练数据L，剩余1/4作为验证数据Y，使用故障数据L对BP神经网络和支持向量回归两种算法进行训练，得到预测训练模型h1和h2；
S412、从未标记数据库中选取未标记数据yi：，使用第j种算法对其进行预测，得到预测输出值，并将第i组未标记数据及其预测的输出与故障数据L组合成新的故障数据；
S413、使用步骤S412得到的新的故障数据对模型重新训练得到新的网络hj’；
S414、使用验证数据Y进行验证，将Y输入hj，计算输出值与实际值之间的均方差，记为eij，再将Y输入hj’，计算输出值与实际值之间的均方差，记为ej’，然后计算eij-eij’的值，记为Eij；
S415、迭代执行步骤S412-S414直到数据库中没有未标记数据为止；
S416、判断所有的Eij中是否有大于0的值，如没有则重新回到S411步对另一种算法进行训练；若存在大于0的值则选择最大的Eij对应的那组未标记数据和网络的输出数据，将其作为伪标记数据，与原有的故障数据组合成新的故障数据组Lj’，并将第i组未标记数据从未标记数据库中移除；对另一种算法进行S412-S415的训练，进入下一步；
S417、判断两种算法训练后L值是否都有更新，若是则交换两种算法的故障数据作为彼此的故障数据，重新进行步骤S411-S416，直到达到迭代次数T为止，若否则结束训练，进行下一步；
S418、对得到的两种网络进行权值优化，两种算法输出值加权后作为最后的网络输出，得到轴承剩余寿命预测值。
本发明的有益效果如下：
本发明所述技术方案充分的考虑到不同故障下轴承衰退模式之间的不同，提高了剩余寿命预测的精度，建立了一整套的轴承故障状态检测与预测体系，完善了对的轴承健康状态监测。使用循环平稳及多变量预测模型进行故障诊断，充分的利用了轴承振动信号的二阶循环特性及特征值之间的内在关系。使用基于半监督的协同训练解决了故障数据难以获取的问题。为轴承后期的维修提供了支持。
附图说明
下面结合附图对本发明的具体实施方式作进一步详细的说明；
图1示出基于故障诊断的轨道交通车辆轴承剩余寿命预测方法流程图。
图2示出轴承振动信号的时域波形图。
图3示出轴承振动信号的频域波形图。
图4示出轴承振动信号的能量波形图。
图5示出轴承振动信号的统计特征波形图。
图6示出特征值融合后的主元成分示意图。
图7示出部分特征值波动点标记示意图。
图8示出不同特征频率下的谱相关密度函数最大值示意图。
图9示出不同特征频率下的谱相关密度函数切片示意图。
图10示出基于BP神经网络、SVR及半监督算法的剩余寿命预测曲线图。
图11示出不同未标记数据下半监督算法的剩余寿命预测曲线图。
图12示出故障后基于不同未标记数据下的半监督算法的剩余寿命预测曲线放大图。
具体实施方式
为了更清楚地说明本发明，下面结合优选实施例和附图对本发明做进一步的说明。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解，下面所具体描述的内容是说明性的而非限制性的，不应以此限制本发明的保护范围。
本发明公开的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法需要对实时采集的振动数据进行分析，通过多角度特征综合处理，提取出能明显表现故障的特征。利用旋转机械振动信号具有循环平稳特性进行故障检测，根据不同故障状态下特征值之间的内在关系使用多变量模型进行模式识别判断故障类型，最后将故障诊断的结果输入不同剩余寿命预测模型，剩余寿命预测主要使用半监督协同训练的算法。
本发明公开的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法相对于已有的剩余寿命预测算法，充分的考虑到了故障类型对于寿命预测结果的影响，结合故障诊断结果对轴承进行剩余寿命预测。对于传统的时域、频域以及时频域故障检测方法中，忽略了旋转机械具有不同于其他机械的运转模式，而这种循环往复的运行方式将产生具周期性的振动信号，包含了大量的状态信息。轴承振动信号的二阶统计量具有周期性变化的特性，因此基于轴承的循环平稳特性，结合经典周期图法的计算方法，使用基于周期图估计循环谱相关密度函数进行轴承的早期故障检测。在故障分离方面充分考虑到不同故障状态下特征值之间具有不同内在关系这一特性，使用多变量预测模型的模式识别法，解决了现有的基于神经网络和支持向量机法在参数设置和核函数选取方面的限制。
本发明公开的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，考虑到目前全生命周期数据难以获取，而大量的未故障状态下数据较容易采集的情况，充分利用未故障状态下的数据信息，使用半监督原理提高算法的预测能力。此外基于现有算法之间的差异，使用BP神经网络及支持向量回归两种算法进行协同训练，防止了算法训练过程中的误差传递。剩余寿命算法中的输入为特征值，输出为剩余寿命值，故障数据已知其输出，记为标记数据；未故障状态下数据的输出未知，记为未标记数据(UL)。
本发明公开的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法以试验平台下的滚动轴承垂向振动加速度为分析对象，主要研究轴承基于故障诊断下的剩余寿命预测。前期使用特征值多角度提取与融合技术，反映轴承的状态特征；通过对循环谱相关密度函数的计算，进行故障的早期检测，结合故障状态下的数据对故障轴承进行模式识别。最后将故障诊断的结果输入剩余寿命预测算法中，对故障早期时刻某种故障模式下的轴承进行剩余寿命预测，实现了轴承可使用时间的预测。
本发明公开的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法包括如下步骤：
S100、进行多角度特征值提取与融合：
在特征提取方面先使用局部均值分解(local mean decomposition,LMD)法将轴承振动信号分解为若干PF单分量；对分离出的单分量进行时域、频域、能量及统计量的多角度特征值提取；最后将提取到的多角度特征值使用主成分分析法(Principal ComponentAnalysis,PCA)降维与融合。
步骤S100的具体过程为：
首先，进行局部均值分解：
S111、寻找轴承振动信号(原始信号)x(t)所有局部极值点ni，将上下极值点分别用三次样条曲线进行连接，得到信号的上、下包络线Emax、Emin，计算局部均值函数m11(t)和包络函数c11(t)：
S112、从x(t)中分离出局部均值函数m11(t)，得到：
h11(t)＝x(t)-m11(t)；
S113、对h11(t)解调，将h11(t)除以包络函数c11(t)得到：
S114、按照S111中的方法求出s11(t)所对应的包络函数c12(t)，若包络函数c12(t)＝1，则判断s11(t)为纯调频函数；若包络函数c12(t)≠1，则迭代执行上述步骤S111-S113，直至s1n(t)的包络估计函数c1(n+1)(t)＝1，则有：
S115、将进行迭代处理过程中产生的所有包络函数相乘得包络信号：
S116、令包络信号c1(t)与纯调频函数s1n(t)相乘，得到x(t)的第一个PF分量PF1：
PF1(t)＝c1(t)·s1n(t)；
S117、从x(t)中分离出第一个分量PF1后得到一个新的待分解信号r1(t)，将新的待分解信号替代轴承振动信号迭代执行步骤S111-S116，直至新的待分解信号rk(t)为一单调函数为止，k代表迭代次数，即
至此，x(t)被分解为k个PF分量及一个单调函数rk(t)之和：
在进行局部均值分解之后，进行多角度特征值提取：
在多角度特征提取中本发明提取18个特征值，包括：
时域特征值包括反映运行状态全局特性并判断轴承是否故障的均方值、方差、峰值、峭度、偏度、脉冲因子、峰值因子和裕度因子等；
频域特征值包括描述信号频域特征的变化并反映轴承的故障部位的均方频率、频谱重心、频率方差和频率标准差等；
统计特征参数特征值包括威布尔形状参数、威布尔尺度参数、伽马形状参数和伽马尺度参数等；
能量特征值包括香农熵(Shannon entropy)和能量矩
在进行多角度特征值提取之后，进行多角度特征值降维与融合：
特征值降维与融合方面使用PCA法，主要过程如下：
S131、设多角度特征值构成的数据的样本数为N，每个样本包含P个特征向量，样本矩阵为XNxP，对数据进行标准化，以降低不同特征值不在同一数量级的问题：
其中，i＝1,2,…,n；；j＝1,2,…,p；；
S132、计算特征向量的协方差矩阵R，并计算其特征值和特征向量：
R＝YTY/(N-1)
其中，Y是标准化后的数据，计算出R的特征值分别为λ1≥λ2≥…≥λp，以及其对应的特征向量为αi＝(αi1,αi2,…,αip)T,i＝1,2,…p；
S133、设正交空间中前k个主元分量为y1,y2…yk，计算累计方差贡献率：
至此，完成了将P维数据降至K维数据，实现了多角度特征值降维与融合。
S200、进行基于循环平稳理论的故障诊断，具体过程如下：
S211、对于为非平稳随机信号的轴承振动信号x(t)，x(t)的时变自相关函数的表示为：
其中，τ为时延因子，E{·}表示统计平均，*表示复共轭；
S212、设Rx(t,τ)以T为周期，用样本平均代替统计平均，将时变自相关函数表示为：
将上式用傅里叶级数展开表示为：
其中α＝m/T(m∈Z)为循环频率，其傅里叶系数为
S213、设T0＝(2N+1)T，通过整理可以得到循环自相关函数为：
R(τ,α)＝<x(t+τ/2)x*(t-τ/2)e-j2παt>t；
S214、得到循环谱相关密度函数并根据循环谱相关密度函数进行基于循环平稳理论的故障诊断，循环谱相关密度函数(Spectral CorrelationDensity，简称SCD)是自相关函数关于时延的Fourier变换，如下式所示：
其中，f为谱频率。
S300、通过进行多变量模型训练和故障类型预测实现基于多变量模型的故障分离：
对于VPMCD而言可以选择的模型如下，
线性模型(L)：
线性交互模型(LI)
二次交互模型(QI)
二次模型(Q)
VPMCD具体操作步可分为进行多变量模型训练和故障类型预测两步，模型训练阶段需要通过对比分析，从可选择的模型中选择最优的变量的预测模型。
步骤S300的具体过程如下：
首先进行多变量模型训练：
S311、针对具有H种故障类型的数据，样本总数为N，每种故障下样本数为[n1,n2,…nH]；
S312、对于不同故障状态的所有样本进行振动数据特征提取，使用拉普拉斯分值法对所有特征值进行选取，每种故障都选出p个最优特征量，分别为X＝[X1,X2,…,Xp]；
S313、设模型类型m＝1，0＜m＜4，模型阶数r＝1，0＜r＜p-1，设h＝1；
S314、对第h(1≤h≤H)类训练样本执行：选择特征量Xi为被预测变量，选择r个特征量Xj(j≠i)为预测变量，则共有种选择方案，即对于被预测变量来说可建立中预测模型；
S315、第h类故障类型共有nh个训练样本，对于每种组合方式，被预测变量Xi都可以建立nh个方程，使用最小二乘法对模型bo,bj,bjj,bjk进行参数估计，得到每种组合下的被预测变量的数学方程
S316、计算每种组合下所有样本对Xi的估计误差和其中l表示第l个训练样本，共有个SSE值，选择最小的SSE对应的预测模型为该种故障，则模型类型为m且阶数为r的条件下，变量Xi的预测模型已确定；
S317、令h＝h+1，迭代执行步骤S313-S316，直到h＝H；
S318、令r＝r+1，迭代执行S314-S317直到r＝P-1；
S319、令m＝m+1，迭代执行S314-S318直到m＝4，得到了各种模型类型及各种阶次条件下的预测模型即每个m和r的条件下都会得到H*P个组成矩阵，然后将训练变量回代各个方程，选择最小误差值下的m和r对应的矩阵为本次的训练输出模型。
之后，进行故障类型预测：
S321、对预测样本计算同样的p个最优特征量，分别为X＝[X1,X2,…,Xp]；
S322、将特征值带入训练好的数学模型，得到预测的估计变量值其中h＝[1,2,…,H]，i＝[1,2,…,p]；
S323、计算每种类型下，所有估计变量的误差值之和SSEh最小误差值的状态类型为该待预测数据的故障类型。
S400、基于半监督算法，协同使用BP神经网络和支持向量回归两种算法进行基于故障诊断的剩余寿命预测：
使用BP神经网络和支持向量回归两种算法进行半监督的协同循环训练，主要步骤如下：将故障数据等间隔取值，3/4作为训练数据L，剩余1/4作为验证数据Y。步骤S400的具体过程为：
S411、初始化：使用故障数据L对BP神经网络和支持向量回归两种算法进行训练，得到预测训练模型h1和h2；
S412、从未标记数据库中选取未标记数据yi：，使用第j种算法对其进行预测，得到预测输出值，并将第i组未标记数据及其预测的输出与故障数据L组合成新的故障数据，其中第j种算法为神经网络算法或支持向量回归算法，用第j种算法来表示是为了后面参数表述方便；
S413、使用步骤S412得到的新的故障数据对模型重新训练得到新的网络hj’；
S414、使用验证数据Y进行验证，将Y输入hj，计算输出值与实际值之间的均方差，记为eij，再将Y输入hj’，计算输出值与实际值之间的均方差，记为ej’，然后计算eij-eij’的值，记为Eij；
S415、迭代执行步骤S412-S414直到数据库中没有未标记数据为止；
S416、判断所有的Eij中是否有大于0的值，如没有则重新回到S411步对另一种算法进行训练；若存在大于0的值则选择最大的Eij对应的那组未标记数据和网络的输出数据，将其作为伪标记数据，与原有的故障数据组合成新的故障数据组Lj’，并将第i组未标记数据从未标记数据库中移除；对另一种算法进行S412-S415的训练，进入下一步；
S417、判断两种算法训练后L值是否都有更新，若是则交换两种算法的故障数据作为彼此的故障数据，重新进行步骤S411-S416，直到达到迭代次数T为止，若否则结束训练，进行下一步；
S418、对得到的两种网络进行权值优化，两种算法输出值加权后作为最后的网络输出，得到轴承剩余寿命预测值。
下面代入具体的轴承振动信号对本发明公开的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法作进一步地说明。
本发明公开的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法包括如下步骤：
S100、进行多角度特征值提取与融合：
采集轴承全生命周期数据，使用LMD法对轴承振动信号进行分解，对得到的单分量进行多角度特征值提取，以PF1分量为例，如图2所示为轴承时域特征值、图3所示为轴承频域特征值、图4所示为轴承能量特征值、图5所示为轴承统计特征值，共提取18种特征值。通过全生命周期的特征值图可以看到每个特征值都能够反映出轴承的运行状态，即在轴承运行初期特征值均显示平稳状态，没有较大波动，但当轴承运行一段时间后特征值开始有波动，呈现向上(向下)的趋势且斜率逐渐增大(减小)
对计算的特征值进行PCA融合降维，取贡献率为95％的主元成分，结果如图6所示，共得到六维主元。可以看出这些分量都完好的保存了轴承数据退化的状态。
S200、进行基于循环平稳理论的故障诊断：
以均方值、方差、峰值和峭度四个时域特征指标为例，标记出有微弱波动的时刻如图7所示。四个特征值的时刻各不相同，最早的在99h处，最晚的在109h处，相差间隔为10小时，对于城轨列车轴承来说，10小时内一直处于故障状态运行，会带来极大的危险。此外四个特征值在118h处均出现较大波动，说明轴承已经处于较为严重的故障状态。通过计算，轴承外圈故障频率为236.4Hz，内圈故障频率为296.8Hz，转频为33.33Hz，滚珠故障频率为139.9Hz，保持架故障频率为14.8Hz。计算各特征频率处的最大谱相关密度函数值如图8所示，从图中可以看出，只有外圈故障频率处在89.67处有明显的波动，其他频率处一直保持平稳状态。
使用89.67处轴承振动信号进行谱切片分析，如图9所示。外圈故障频率附近谱密度相关切片图中237.4Hz处的幅值最大，此处237.4Hz与理论的外圈故障频率236.4Hz存在一定的差异，主要是由于转速的波动和实际承载间的变化造成实际故障频率与理论故障频率的差异；同理将内圈、滚动体、保持架故障频率及转频附近的谱密度相关切片中能量最大的切片取出。将取出的最大能量切片组合，图中可以看出只有外圈故障频率处的切片幅值最高，因此判定该轴承外圈发生故障。
S300、通过进行多变量模型训练和故障类型预测实现基于多变量模型的故障分离：
选择特波形因子，峭度指标，功率谱重心指标，均值频率，均方谱作为特征值。其中训练数据每种状态40个样本，测试数据每种状态10个样本。数据准备好后根据VPMCD的流程进行模型训练。再将89.67处的振动信号同样进行特征值计算，输入训练好的数学模型中，计算每种状态下5个特征值的估计误差总和，对每种状态下误差之和进行对比，选择误差最小的一种状态作为该测试数据的真正状态，测试结果，即轴承89.67小时处的信号故障分离结果如表1所示，结果与实际相符。对10个测试样本进行预测，使用BP神经网络和支持向量机算法与VPMCD在训练时间、测试结果的准确度上进行对比分析如表2所示。
表1
表2
S400、协同使用BP神经网络和支持向量回归两种算法进行基于故障诊断的剩余寿命预测：
使用经过粒子群优化后的BP神经网络和支持向量回归作为协同训练的两种算法，使用另外8组未标记数据以及一组外圈故障数据最为算法的训练样本，使用检测到故障的数据作为测试样本。
使用两点的特征值作为网络模型的输入，两点即为当前点对应的特征值和前一时刻监测点对应的特征值，维数为12；输出为剩余寿命值，维数为1。
对单独使用BP神经网络和支持向量回归时的预测结果与使用8个未标记数据时的预测结果进行对比，如图10所示；此外针对不同的未标记数据个数对于算法精度的影响进行对比如图11所示，图12所示是将图11在故障后(85h之后)的预测图形放大图。使用均方根误差、平均绝对误差、希尔不等系数和平均相对变动值作为网络的评价指标，计算结果如表3所示。通过对比可以说明使用未标记数据对剩余寿命预测的算法精度有一定的提升，随着未标记数据的增多，预测精度有所升高，且波动情况相对平缓。
表3
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (7)

1.一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，其特征在于，该方法包括如下步骤：
S100、使用局部均值分解法将轴承振动信号分解为若干PF单分量，对分离出的单分量进行时域、频域、能量及统计量的多角度特征值提取，将提取到的多角度特征值使用主成分分析法降维与融合；
S200、进行基于循环平稳理论的故障诊断；
S300、通过进行多变量模型训练和实现基于多变量模型的故障分离；
S400、基于半监督算法，协同使用BP神经网络和支持向量回归两种算法进行基于故障诊断的剩余寿命预测。
2.根据权利要求1所述的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，其特征在于，所述使用局部均值分解法将轴承振动信号分解为若干PF单分量进一步包括：
S111、寻找轴承振动信号x(t)所有局部极值点ni，将上下极值点分别用三次样条曲线进行连接，得到信号的上、下包络线Emax、Emin，计算局部均值函数m11(t)和包络函数c11(t)：
<mrow> <msub> <mi>m</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>E</mi> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>E</mi> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mn>2</mn> </mfrac> </mrow>
<mrow> <msub> <mi>c</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <msub> <mi>E</mi> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>E</mi> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mn>2</mn> </mfrac> </mrow>
S112、从x(t)中分离出局部均值函数m11(t)，得到：
h11(t)＝x(t)-m11(t)；
S113、对h11(t)解调，将h11(t)除以包络函数c11(t)得到：
<mrow> <msub> <mi>s</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>h</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>c</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>;</mo> </mrow>
S114、按照S111中的方法求出s11(t)所对应的包络函数c12(t)，若包络函数c12(t)＝1，则判断s11(t)为纯调频函数；若包络函数c12(t)≠1，则迭代执行上述步骤S111-S113，直至s1n(t)的包络估计函数c1(n+1)(t)＝1，则有：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>h</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>m</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>h</mi> <mn>12</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>s</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>m</mi> <mn>12</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>h</mi> <mrow> <mn>1</mn> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>s</mi> <mrow> <mn>1</mn> <mrow> <mo>(</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>m</mi> <mrow> <mn>1</mn> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>s</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>h</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>c</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>s</mi> <mn>12</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>h</mi> <mn>12</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>c</mi> <mn>12</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>s</mi> <mrow> <mn>1</mn> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>h</mi> <mrow> <mn>1</mn> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>c</mi> <mrow> <mn>1</mn> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> </mrow>
S115、将进行迭代处理过程中产生的所有包络函数相乘得包络信号：
<mrow> <msub> <mi>c</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>c</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <msub> <mi>c</mi> <mn>12</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>...</mo> <msub> <mi>c</mi> <mrow> <mn>1</mn> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Pi;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>c</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
S116、令包络信号c1(t)与纯调频函数s1n(t)相乘，得到x(t)的第一个PF分量PF1：
PF1(t)＝c1(t)·s1n(t)；
S117、从x(t)中分离出第一个分量PF1后得到一个新的待分解信号r1(t)，将新的待分解信号替代轴承振动信号迭代执行步骤S111-S116，直至新的待分解信号rk(t)为一单调函数为止，k代表迭代次数，即
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>r</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>PF</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>r</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>r</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>PF</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>r</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>PF</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> </mrow>
至此，x(t)被分解为k个PF分量及一个单调函数rk(t)之和：
<mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>v</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>k</mi> </munderover> <msub> <mi>PF</mi> <mi>v</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>r</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>.</mo> </mrow>
3.根据权利要求1所述的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，其特征在于，所述多角度特征值包括时域特征值、频域特征值、统计特征参数特征值和能量特征值；所述时域特征值包括均方值、方差、峰值、峭度、偏度、脉冲因子、峰值因子和裕度因子；所述频域特征值包括均方频率、频谱重心、频率方差和频率标准差；所述统计特征参数特征值包括威布尔形状参数、威布尔尺度参数、伽马形状参数和伽马尺度参数；所述能量特征值包括香农熵和能量矩。
4.根据权利要求1所述的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，其特征在于，所述将提取到的多角度特征值使用主成分分析法降维与融合进一步包括：
S131、设多角度特征值构成的数据的样本数为N，每个样本包含P个特征向量，样本矩阵为XNxP，对数据进行标准化：
<mrow> <msub> <mi>k</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>x</mi> <mo>&amp;OverBar;</mo> </mover> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>&amp;sigma;</mi> <mi>j</mi> </msub> </mrow>
<mrow> <msub> <mi>&amp;sigma;</mi> <mi>j</mi> </msub> <mo>=</mo> <msqrt> <mrow> <mn>1</mn> <mo>/</mo> <mrow> <mo>(</mo> <mi>N</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msqrt> <mo>&amp;CenterDot;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>x</mi> <mo>&amp;OverBar;</mo> </mover> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow>
其中，i＝1,2,…,n；；j＝1,2,…,p；；
S132、计算特征向量的协方差矩阵R，并计算其特征值和特征向量：
R＝YTY/(N-1)
其中，Y是标准化后的数据，计算出R的特征值分别为λ1≥λ2≥…≥λp，以及其对应的特征向量为αi＝(αi1,αi2,…,αip)T,i＝1,2,…p；
S133、设正交空间中前k个主元分量为y1,y2…yk，计算累计方差贡献率：
<mrow> <mi>&amp;theta;</mi> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>k</mi> </munderover> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mo>/</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msub> <mi>&amp;lambda;</mi> <mi>j</mi> </msub> </mrow>
至此，完成了将P维数据降至K维数据，实现了多角度特征值降维与融合。
5.根据权利要求1所述的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，其特征在于，步骤S200进一步包括：
S211、对于为非平稳随机信号的轴承振动信号x(t)，x(t)的时变自相关函数的表示为：
<mrow> <msub> <mi>R</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>E</mi> <mo>{</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>+</mo> <mfrac> <mi>&amp;tau;</mi> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <msup> <mi>x</mi> <mo>*</mo> </msup> <mrow> <mo>(</mo> <mi>t</mi> <mo>-</mo> <mfrac> <mi>&amp;tau;</mi> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <mo>}</mo> </mrow>
其中，τ为时延因子，E{·}表示统计平均，*表示复共轭；
S212、设Rx(t,τ)以T为周期，用样本平均代替统计平均，将时变自相关函数表示为：
<mrow> <munder> <mi>lim</mi> <mrow> <mi>N</mi> <mo>&amp;RightArrow;</mo> <mi>&amp;infin;</mi> </mrow> </munder> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>N</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mo>-</mo> <mi>N</mi> </mrow> <mi>N</mi> </munderover> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>+</mo> <mi>n</mi> <mi>T</mi> <mo>+</mo> <mfrac> <mi>&amp;tau;</mi> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <msup> <mi>x</mi> <mo>*</mo> </msup> <mrow> <mo>(</mo> <mi>t</mi> <mo>+</mo> <mi>n</mi> <mi>T</mi> <mo>-</mo> <mfrac> <mi>&amp;tau;</mi> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow>
将上式用傅里叶级数展开表示为：
<mrow> <msub> <mi>R</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>&amp;infin;</mi> </mrow> <mi>&amp;infin;</mi> </munderover> <msubsup> <mi>R</mi> <mi>x</mi> <mi>&amp;alpha;</mi> </msubsup> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>&amp;alpha;</mi> <mi>t</mi> </mrow> </msup> </mrow>
其中α＝m/T(m∈Z)为循环频率，其傅里叶系数为
<mrow> <msubsup> <mi>R</mi> <mi>x</mi> <mi>&amp;alpha;</mi> </msubsup> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mi>T</mi> </mfrac> <msubsup> <mo>&amp;Integral;</mo> <mrow> <mo>-</mo> <mi>T</mi> <mo>/</mo> <mn>2</mn> </mrow> <mrow> <mi>T</mi> <mo>/</mo> <mn>2</mn> </mrow> </msubsup> <msub> <mi>R</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>&amp;alpha;</mi> <mi>t</mi> </mrow> </msup> <mi>d</mi> <mi>t</mi> <mo>;</mo> </mrow>
S213、设T0＝(2N+1)T，通过整理可以得到循环自相关函数为：
R(τ,α)＝<x(t+τ/2)x*(t-τ/2)e-j2παt>t；
S214、得到循环谱相关密度函数并根据循环谱相关密度函数进行基于循环平稳理论的故障诊断，循环谱相关密度函数为：
<mrow> <mi>S</mi> <mrow> <mo>(</mo> <mi>f</mi> <mo>,</mo> <mi>&amp;alpha;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&amp;Integral;</mo> <mrow> <mo>-</mo> <mi>&amp;infin;</mi> </mrow> <mi>&amp;infin;</mi> </msubsup> <msub> <mi>R</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>&amp;alpha;</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>f</mi> <mi>&amp;tau;</mi> </mrow> </msup> <mi>d</mi> <mi>&amp;tau;</mi> </mrow>
其中，f为谱频率。
6.根据权利要求1所述的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，其特征在于，步骤S300进一步包括：
进行多变量模型训练：
S311、针对具有H种故障类型的数据，样本总数为N，每种故障下样本数为[n1,n2,…nH]；
S312、对于不同故障状态的所有样本进行振动数据特征提取，使用拉普拉斯分值法对所有特征值进行选取，每种故障都选出p个最优特征量，分别为X＝[X1,X2,…,Xp]
S313、设模型类型m＝1，0＜m＜4，模型阶数r＝1，0＜r＜p-1，设h＝1；
S314、对第h(1≤h≤H)类训练样本执行：选择特征量Xi为被预测变量，选择r个特征量Xj(j≠i)为预测变量，则共有种选择方案，即对于被预测变量来说可建立中预测模型；
S315、第h类故障类型共有nh个训练样本，对于每种组合方式，被预测变量Xi都可以建立nh个方程，使用最小二乘法对模型bo,bj,bjj,bjk进行参数估计，得到每种组合下的被预测变量的数学方程
S316、计算每种组合下所有样本对Xi的估计误差和其中l表示第l个训练样本，共有个SSE值，选择最小的SSE对应的预测模型为该种故障，则模型类型为m且阶数为r的条件下，变量Xi的预测模型已确定；
S317、令h＝h+1，迭代执行步骤S313-S316，直到h＝H；
S318、令r＝r+1，迭代执行S314-S317直到r＝P-1；
S319、令m＝m+1，迭代执行S314-S318直到m＝4，得到了各种模型类型及各种阶次条件下的预测模型即每个m和r的条件下都会得到H*P个组成矩阵，然后将训练变量回代各个方程，选择最小误差值下的m和r对应的矩阵为本次的训练输出模型；
进行故障类型预测：
S321、对预测样本计算同样的p个最优特征量，分别为X＝[X1,X2,…,Xp]；
S322、将特征值带入训练好的数学模型，得到预测的估计变量值其中h＝[1,2,…,H]，i＝[1,2,…,p]；
S323、计算每种类型下，所有估计变量的误差值之和SSEh最小误差值的状态类型为该待预测数据的故障类型。
7.根据权利要求1所述的基于故障诊断的轨道交通车辆轴承剩余寿命预测方法，其特征在于，步骤S400进一步包括：
S411、将故障数据等间隔取值，3/4作为训练数据L，剩余1/4作为验证数据Y，使用故障数据L对BP神经网络和支持向量回归两种算法进行训练，得到预测训练模型h1和h2；
S412、从未标记数据库中选取未标记数据yi：，使用第j种算法对其进行预测，得到预测输出值，并将第i组未标记数据及其预测的输出与故障数据L组合成新的故障数据；
S413、使用步骤S412得到的新的故障数据对模型重新训练得到新的网络hj’；
S414、使用验证数据Y进行验证，将Y输入hj，计算输出值与实际值之间的均方差，记为eij，再将Y输入hj’，计算输出值与实际值之间的均方差，记为ej’，然后计算eij-eij’的值，记为Eij；
S415、迭代执行步骤S412-S414直到数据库中没有未标记数据为止；
S416、判断所有的Eij中是否有大于0的值，如没有则重新回到S411步对另一种算法进行训练；若存在大于0的值则选择最大的Eij对应的那组未标记数据和网络的输出数据，将其作为伪标记数据，与原有的故障数据组合成新的故障数据组Lj’，并将第i组未标记数据从未标记数据库中移除；对另一种算法进行S412-S415的训练，进入下一步；
S417、判断两种算法训练后L值是否都有更新，若是则交换两种算法的故障数据作为彼此的故障数据，重新进行步骤S411-S416，直到达到迭代次数T为止，若否则结束训练，进行下一步；
S418、对得到的两种网络进行权值优化，两种算法输出值加权后作为最后的网络输出，得到轴承剩余寿命预测值。
CN201710160858.8A 2017-03-17 2017-03-17 一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法 Active CN107144428B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710160858.8A CN107144428B (zh) 2017-03-17 2017-03-17 一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710160858.8A CN107144428B (zh) 2017-03-17 2017-03-17 一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法

Publications (2)

Publication Number Publication Date
CN107144428A true CN107144428A (zh) 2017-09-08
CN107144428B CN107144428B (zh) 2019-02-22

Family

ID=59783427

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710160858.8A Active CN107144428B (zh) 2017-03-17 2017-03-17 一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法

Country Status (1)

Country Link
CN (1) CN107144428B (zh)

Cited By (30)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107831013A (zh) * 2017-10-11 2018-03-23 温州大学 一种利用概率主分量分析增强循环双谱的轴承故障诊断方法
CN108304960A (zh) * 2017-12-29 2018-07-20 中车工业研究院有限公司 一种轨道交通设备故障诊断方法
CN108776818A (zh) * 2018-06-05 2018-11-09 北京航空航天大学 轴承特征提取方法、轴承剩余寿命预测方法和装置
CN108896312A (zh) * 2018-08-08 2018-11-27 国电联合动力技术有限公司 一种风电主轴承故障预测和寿命评估系统及方法
CN108985279A (zh) * 2018-08-28 2018-12-11 上海仁童电子科技有限公司 多功能车辆总线mvb波形的故障诊断方法及装置
CN109142428A (zh) * 2018-04-13 2019-01-04 权奥博 一种动车轨道强度智能检测系统
CN109190830A (zh) * 2018-09-11 2019-01-11 四川大学 基于经验分解和组合预测的能源需求预测方法
CN109242215A (zh) * 2018-10-26 2019-01-18 浙江工业大学之江学院 一种基于粒子群-支持向量机的旋转机械设备运行工况预测方法
CN109472241A (zh) * 2018-11-14 2019-03-15 上海交通大学 基于支持向量回归的燃机轴承剩余使用寿命预测方法
CN109580224A (zh) * 2018-12-28 2019-04-05 北京中科东韧科技有限责任公司 滚动轴承故障实时监测方法
CN109614850A (zh) * 2018-10-26 2019-04-12 武汉长天铁路技术有限公司 基于L-M的More算法的轨道不平顺谱拟合方法
CN109636250A (zh) * 2019-01-17 2019-04-16 长安大学 一种危险货物卡车生存概率和危险概率的预测方法
CN110320018A (zh) * 2019-07-12 2019-10-11 北京交通大学 一种基于二阶循环平稳特性的旋转机械复合故障诊断方法
CN110333265A (zh) * 2019-07-11 2019-10-15 中国中车股份有限公司 一种预测机车发动机排气波纹管剩余寿命的方法及系统
CN110608885A (zh) * 2019-09-09 2019-12-24 天津工业大学 一种用于滚动轴承内圈磨损故障诊断及趋势预测的方法
CN110688959A (zh) * 2019-09-27 2020-01-14 上海特金信息科技有限公司 无人机信号识别方法、装置、电子设备与存储介质
CN110766264A (zh) * 2019-07-31 2020-02-07 贵州电网有限责任公司 一种传输电缆健康状态监测评估方法
CN110987396A (zh) * 2019-12-13 2020-04-10 三一重型装备有限公司 一种用于采煤机摇臂的智能故障诊断及寿命预测方法
CN111194431A (zh) * 2018-09-14 2020-05-22 西安大医集团有限公司 放疗设备的状态诊断方法、装置、系统及存储介质
CN111220388A (zh) * 2020-02-04 2020-06-02 珠海市华星装备信息科技有限公司 一种基于时域分析的滚动轴承故障诊断方法
CN111222203A (zh) * 2018-11-08 2020-06-02 上海仪电（集团）有限公司中央研究院 一种轴承使用寿命模型创建及其预测方法
CN111259737A (zh) * 2020-01-08 2020-06-09 科大讯飞股份有限公司 车辆方向盘故障的预测方法、装置、电子设备及存储介质
US20210056780A1 (en) * 2019-08-22 2021-02-25 GM Global Technology Operations LLC Adaptive fault diagnostic system for motor vehicles
CN112580153A (zh) * 2020-12-29 2021-03-30 成都运达科技股份有限公司 一种车辆走行部监测部件健康状态管理系统及方法
CN112669262A (zh) * 2020-12-08 2021-04-16 上海交通大学 一种电机轮轴震动异常检测与预测系统与方法
CN112884210A (zh) * 2021-01-31 2021-06-01 中国人民解放军63963部队 一种基于模糊聚类的车辆健康管理系统架构优化方法
WO2021117752A1 (ja) * 2019-12-11 2021-06-17 Ntn株式会社 転がり軸受の状態監視方法及び転がり軸受の状態監視装置
CN113125155A (zh) * 2021-04-29 2021-07-16 浙江陀曼云计算有限公司 工业轴承检测方法及系统
CN113776834A (zh) * 2021-10-11 2021-12-10 山东大学 基于离散余弦循环谱相干的滚动轴承故障诊断方法及系统
CN114757239A (zh) * 2022-06-15 2022-07-15 浙江大学 基于数据增强和胶囊神经网络的风机故障可迁移诊断方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2177891A2 (de) * 2008-10-13 2010-04-21 Chemserv Industrie Service GmbH Vorrichtung zum Überwachen der Wälzlager von rotierenden Maschinen einer Produktionsanlage
CN102435436A (zh) * 2011-11-24 2012-05-02 电子科技大学 风扇轴承状态退化评估方法
CN102721545A (zh) * 2012-05-25 2012-10-10 北京交通大学 一种基于多特征参量的滚动轴承故障诊断方法
KR101492090B1 (ko) * 2013-11-21 2015-02-10 이선휘 구름베어링의 잔여수명 예측방법
CN105938468A (zh) * 2016-06-07 2016-09-14 北京交通大学 一种滚动轴承的故障诊断方法
CN106053066A (zh) * 2016-05-23 2016-10-26 华东交通大学 基于经验模态分解和逻辑回归的滚动轴承性能退化评估方法
CN106248380A (zh) * 2016-09-09 2016-12-21 芜湖能盟信息技术有限公司 一种轴承寿命预测试验方法及其系统

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2177891A2 (de) * 2008-10-13 2010-04-21 Chemserv Industrie Service GmbH Vorrichtung zum Überwachen der Wälzlager von rotierenden Maschinen einer Produktionsanlage
CN102435436A (zh) * 2011-11-24 2012-05-02 电子科技大学 风扇轴承状态退化评估方法
CN102721545A (zh) * 2012-05-25 2012-10-10 北京交通大学 一种基于多特征参量的滚动轴承故障诊断方法
KR101492090B1 (ko) * 2013-11-21 2015-02-10 이선휘 구름베어링의 잔여수명 예측방법
CN106053066A (zh) * 2016-05-23 2016-10-26 华东交通大学 基于经验模态分解和逻辑回归的滚动轴承性能退化评估方法
CN105938468A (zh) * 2016-06-07 2016-09-14 北京交通大学 一种滚动轴承的故障诊断方法
CN106248380A (zh) * 2016-09-09 2016-12-21 芜湖能盟信息技术有限公司 一种轴承寿命预测试验方法及其系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
王玉静: "滚动轴承振动信号特征提取与状态评估方法研究", 《中国博士学位论文全文数据库 工程科技II辑》 *

Cited By (40)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107831013A (zh) * 2017-10-11 2018-03-23 温州大学 一种利用概率主分量分析增强循环双谱的轴承故障诊断方法
CN107831013B (zh) * 2017-10-11 2019-11-22 温州大学 利用概率主分量分析增强循环双谱的轴承故障诊断方法
CN108304960A (zh) * 2017-12-29 2018-07-20 中车工业研究院有限公司 一种轨道交通设备故障诊断方法
CN109142428A (zh) * 2018-04-13 2019-01-04 权奥博 一种动车轨道强度智能检测系统
CN108776818A (zh) * 2018-06-05 2018-11-09 北京航空航天大学 轴承特征提取方法、轴承剩余寿命预测方法和装置
CN108896312A (zh) * 2018-08-08 2018-11-27 国电联合动力技术有限公司 一种风电主轴承故障预测和寿命评估系统及方法
CN108985279A (zh) * 2018-08-28 2018-12-11 上海仁童电子科技有限公司 多功能车辆总线mvb波形的故障诊断方法及装置
CN109190830A (zh) * 2018-09-11 2019-01-11 四川大学 基于经验分解和组合预测的能源需求预测方法
CN109190830B (zh) * 2018-09-11 2021-11-30 四川大学 基于经验分解和组合预测的能源需求预测方法
CN111194431A (zh) * 2018-09-14 2020-05-22 西安大医集团有限公司 放疗设备的状态诊断方法、装置、系统及存储介质
CN109242215A (zh) * 2018-10-26 2019-01-18 浙江工业大学之江学院 一种基于粒子群-支持向量机的旋转机械设备运行工况预测方法
CN109614850A (zh) * 2018-10-26 2019-04-12 武汉长天铁路技术有限公司 基于L-M的More算法的轨道不平顺谱拟合方法
CN111222203A (zh) * 2018-11-08 2020-06-02 上海仪电（集团）有限公司中央研究院 一种轴承使用寿命模型创建及其预测方法
CN109472241A (zh) * 2018-11-14 2019-03-15 上海交通大学 基于支持向量回归的燃机轴承剩余使用寿命预测方法
CN109580224A (zh) * 2018-12-28 2019-04-05 北京中科东韧科技有限责任公司 滚动轴承故障实时监测方法
CN109636250A (zh) * 2019-01-17 2019-04-16 长安大学 一种危险货物卡车生存概率和危险概率的预测方法
CN110333265A (zh) * 2019-07-11 2019-10-15 中国中车股份有限公司 一种预测机车发动机排气波纹管剩余寿命的方法及系统
CN110320018A (zh) * 2019-07-12 2019-10-11 北京交通大学 一种基于二阶循环平稳特性的旋转机械复合故障诊断方法
CN110320018B (zh) * 2019-07-12 2020-08-11 北京交通大学 一种基于二阶循环平稳特性的旋转机械复合故障诊断方法
CN110766264A (zh) * 2019-07-31 2020-02-07 贵州电网有限责任公司 一种传输电缆健康状态监测评估方法
US11551488B2 (en) * 2019-08-22 2023-01-10 GM Global Technology Operations LLC Adaptive fault diagnostic system for motor vehicles
US20210056780A1 (en) * 2019-08-22 2021-02-25 GM Global Technology Operations LLC Adaptive fault diagnostic system for motor vehicles
CN110608885A (zh) * 2019-09-09 2019-12-24 天津工业大学 一种用于滚动轴承内圈磨损故障诊断及趋势预测的方法
CN110688959A (zh) * 2019-09-27 2020-01-14 上海特金信息科技有限公司 无人机信号识别方法、装置、电子设备与存储介质
WO2021117752A1 (ja) * 2019-12-11 2021-06-17 Ntn株式会社 転がり軸受の状態監視方法及び転がり軸受の状態監視装置
CN110987396B (zh) * 2019-12-13 2021-07-30 三一重型装备有限公司 一种用于采煤机摇臂的智能故障诊断及寿命预测方法
CN110987396A (zh) * 2019-12-13 2020-04-10 三一重型装备有限公司 一种用于采煤机摇臂的智能故障诊断及寿命预测方法
CN111259737A (zh) * 2020-01-08 2020-06-09 科大讯飞股份有限公司 车辆方向盘故障的预测方法、装置、电子设备及存储介质
CN111259737B (zh) * 2020-01-08 2023-07-25 科大讯飞股份有限公司 车辆方向盘故障的预测方法、装置、电子设备及存储介质
CN111220388B (zh) * 2020-02-04 2021-11-30 珠海市华星装备信息科技有限公司 一种基于时域分析的滚动轴承故障诊断方法
CN111220388A (zh) * 2020-02-04 2020-06-02 珠海市华星装备信息科技有限公司 一种基于时域分析的滚动轴承故障诊断方法
CN112669262A (zh) * 2020-12-08 2021-04-16 上海交通大学 一种电机轮轴震动异常检测与预测系统与方法
CN112580153B (zh) * 2020-12-29 2022-10-11 成都运达科技股份有限公司 一种车辆走行部监测部件健康状态管理系统及方法
CN112580153A (zh) * 2020-12-29 2021-03-30 成都运达科技股份有限公司 一种车辆走行部监测部件健康状态管理系统及方法
CN112884210A (zh) * 2021-01-31 2021-06-01 中国人民解放军63963部队 一种基于模糊聚类的车辆健康管理系统架构优化方法
CN112884210B (zh) * 2021-01-31 2022-03-15 中国人民解放军63963部队 一种基于模糊聚类的车辆健康管理系统
CN113125155A (zh) * 2021-04-29 2021-07-16 浙江陀曼云计算有限公司 工业轴承检测方法及系统
CN113776834A (zh) * 2021-10-11 2021-12-10 山东大学 基于离散余弦循环谱相干的滚动轴承故障诊断方法及系统
CN114757239B (zh) * 2022-06-15 2022-08-30 浙江大学 基于数据增强和胶囊神经网络的风机故障可迁移诊断方法
CN114757239A (zh) * 2022-06-15 2022-07-15 浙江大学 基于数据增强和胶囊神经网络的风机故障可迁移诊断方法

Also Published As

Publication number Publication date
CN107144428B (zh) 2019-02-22

Similar Documents

Publication Publication Date Title
CN107144428A (zh) 一种基于故障诊断的轨道交通车辆轴承剩余寿命预测方法
KR101316486B1 (ko) 이상 검지 방법 및 시스템
Ragab et al. Remaining useful life prediction using prognostic methodology based on logical analysis of data and Kaplan–Meier estimation
JP5538597B2 (ja) 異常検知方法及び異常検知システム
JP6216242B2 (ja) 異常検知方法およびその装置
Bag et al. An expert system for control chart pattern recognition
CN102778355B (zh) 一种基于emd和pca的滚动轴承状态辨识方法
WO2010095314A1 (ja) 異常検知方法及び異常検知システム
WO2018053935A1 (zh) 一种采用故障模式发生概率的动设备运行状态模糊评价及预测方法
CN111598150B (zh) 一种计及运行状态等级的变压器故障诊断方法
Gandhi et al. Towards data mining based decision support in manufacturing maintenance
Vimal et al. Application of artificial neural network for fuzzy logic based leanness assessment
CN104502103A (zh) 一种基于模糊支持向量机的轴承故障诊断方法
CN103291600B (zh) 一种基于emd-ar和mts的液压泵故障诊断方法
CN106124175A (zh) 一种基于贝叶斯网络的压缩机气阀故障诊断方法
CN111222549A (zh) 一种基于深度神经网络的无人机故障预测方法
Khoshgoftaar et al. Data mining for predictors of software quality
CN109583520B (zh) 一种云模型与遗传算法优化支持向量机的状态评估方法
Yacout Fault detection and diagnosis for condition based maintenance using the logical analysis of data
CN108241901A (zh) 一种基于预测数据的变压器预警评估方法及装置
Faruk et al. Prediction and classification of low birth weight data using machine learning techniques
Garg et al. A two-phase approach for reliability and maintainability analysis of an industrial system
CN106529580A (zh) 结合edsvm的软件缺陷数据关联分类方法
Silva et al. Damage detection in a benchmark structure using AR-ARX models and statistical pattern recognition
CN104832418A (zh) 一种基于局部均值变换和Softmax的液压泵故障诊断方法

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