技术领域
[0001] 本发明涉及故障诊断技术领域，特别是涉及一种基于多维同步最优分解的轴承故障诊断方法。
相关背景技术
[0002] 滚动轴承是旋转机械的核心部件，其工作异常往往会造成设备性能下降乃至故障。因此，滚动轴承故障诊断对于确保机械设备安全运行至关重要。然而，测试信号不仅受到非故障部件和环境的干扰，还受到冗长传递路径的严重衰减。因此，实现轴承故障诊断的关键在于如何从复杂测试信号中提取出与故障相关的周期性冲击。
[0003] 解卷积方法因其能有效抵抗传输路径影响而得到广泛应用。其核心在于设计自适应滤波器，并通过迭代或非迭代方式实现最优滤波器与潜在故障信息的精确匹配。1978年，Wiggins首次提出了最小熵解卷积(Minimum entropy deconvolution，简称MED)，用于从复杂信号中恢复地震脉冲信号的特征。随后，McDonald等人于2012年提出了最大相关峭度解卷积(Maximum correlated kurtosis deconvolution，简称MCKD)，考虑了信号的周期性。2018年，Buzzoni等人提出了最大二阶循环平稳解卷积(Maximum second‑order cyclostationarity blind deconvolution，简称CYCBD)，考虑了旋转机械故障的循环平稳性质。尽管上述方法在大多数情况下简单有效，但它们主要依赖于单通道信息进行故障特征提取，故障信息不足以应对噪声和随机干扰，可能导致漏报。
[0004] 相比之下，分解方法在多组件分离和提取方面更具优势。经验模式分解(Empirical mode decomposition，简称EMD)、经验小波变换(Empirical wavelet transform，简称EWT)和变分模式分解(Variational mode decomposition，简称VMD)等方法被相继提出，但这些传统分解方法并不完全适用于轴承故障诊断，并未考虑典型机械故障特征。为了克服传统分解方法的局限性，Miao等人提出了基于故障特征的模式分解(Feature mode decomposition，简称FMD)，采用相关峭度作为目标函数，具有较高鲁棒性。然而，FMD的输出模式数量仍需手动设置。
[0005] 随着传感器技术的发展，多通道信息输入为故障特征信息挖掘提供了解决方案。EMD、EWT和VMD等方法被扩展到多变量版本，但并未摆脱原始分解方法的缺点。考虑到FMD的优秀特性，Yan等人试图将FMD推广到多元特征模式分解中，通过在单个通道上执行FMD，获得多个组件作为后续卷积神经网络的输入，以学习判别特征。上述多元分解方法处理得到的分解模式众多，需要进一步人工观察来识别故障。此外，人们研究了大量基于数据融合手段进行多通道信号处理的策略，其表现出更优异的性能，但在一定程度上依赖于大量标记数据。这些不足限制了传统分解方法在轴承多通道信号处理和多元故障特征提取与诊断中的应用与推广。
具体实施方式
[0057] 以下，参照附图对本发明的实施方式进行说明。
[0058] 图2为本发明基于多维同步最优分解的轴承故障诊断方法流程图，如图2所示，本发明提供了一种基于多维同步最优分解的轴承故障诊断方法，其包括：
[0059] 步骤S1，获取待测轴承的多通道振动信号Xraw＝[x1 x2 … xO]，其中各通道信号为Txo＝[xo(1) xo(2) … xo(N)] ,o∈O，O为通道数，N为信号的长度，上标T为转置；对所述振动信号进行预处理，包括高频采样、截断和去均值，得到预处理信号；
[0060] 步骤S2，设定滤波器长度为L，滑动窗口长度为W，预迭代次数为I、目标故障类别数为Q，总模式数量为K，理论故障周期为TQ；
[0061] 步骤S3，对O个通道信号进行计算，得到O个上包络集合，W为滑动窗口长度；取各组上包络集合的极小值点作为频谱分割集合，得到各通道初始滤波器截止频率集合；基于Hanning窗构建O组通道滤波器组，各组滤波器数量为Ko,o∈O；
[0062] 所述各通道上包络集合的计算公式如下：
[0063]
[0064] 式中：UpEnvo(n)表示第o个通道信号xo计算所得的上包络集合中第n个长度点位置T的参数值，Yo表示xo的傅里叶频谱，Yo＝[Yo(1) Yo(2) … Yo(N)]，其中Yo(n)表示Yo第n个长度点位置的参数值，fft为傅里叶函数；
[0065] 所述各通道截止频率集合的计算公式如下：
[0066]
[0067] 式中：[flower,fupper]o表示第o组初始滤波器的下截止频率和上截止频率，FreBoundo表示第o个通道信号上包络集合的极小值点集合，Bo为第o个通道信号上包络集合的极小值点数量，FreBoundo(b)表示第o个通道信号上包络集合的极小值点集合第b个长度点位置的参数值；
[0068] 步骤S4，令i＝1，j＝1，i∈I，j∈Q，I为预迭代次数，Q为故障类别数；
[0069] 步骤S5，对i进行判断，若i≤I，则执行S6，若i＞I，则执行S8；
[0070] 步骤S6，基于Q个理论故障周期计算相关峭度，对所述各通道信号进行滤波，对于各滤波器，均选取所述广义瑞利熵等式中特征值的最大值所对应的分解模式作为最终分解模式；
[0071] 所述相关峭度计算公式如下：
[0072]
[0073] 式中：uo,k＝[uo,k(1) uo,k(2) … uo,k(N)]T表示第o个通道的第k个分解模式，o∈O,k∈Ko，uo,k(n)表示第o个通道的第k个分解模式中第n个长度点位置的参数值，M为移位阶数， fo,k为第o个通道的第k个滤波器的向量形式，fo,k＝[fo,kT
(1) fo,k(2) … fo,k(L)]，fo,k(l)表示第o个通道的第k个滤波器中第l个长度点位置的参数值，Tj为第j类故障对应的理论故障周期，j∈Q；
[0074]
[0075] 式中：RXWX和RXX均为加权相关矩阵，xo为第o个通道信号，上标H为共轭转置；
[0077] 广义瑞利熵等式如下：
[0078] fo,kRXWX＝λfo,kRXX；
[0079] 式中：λ为特征值；
[0080] 步骤S7，对所有所述滤波器的系数进行调整，令i＝i+1并返回至S5；
[0081] 步骤S8，对j进行判断，若j≤Q，则执行S9，若j＞Q，则执行S11；
[0082] 步骤S9，对于各通道所述Ko个最终分解模式，均计算所述最终分解模式的相关峭度，得到Ko个相关峭度；计算保留数量P，选取相关峭度的前P个最大值对应的最终分解模式作为所述最优分解模式；
[0083] 所述保留数量的计算公式如下：
[0084] P＝floor(ΣKo/Q·O)；
[0085] 式中：floor为向下取整操作，Ko为各通道的滤波器数量，Q为故障类别数，O为总通道数；
[0086] 步骤S10，令j＝j+1并返回至S8；
[0087] 步骤S11，计算K个所述最优分解模式两两之间的时域相关系数和频域相关系数，选取二者皆大于0.8的所有最优分解模式组合作为筛选分解模式组合；
[0088] 所述时域相关系数的计算公式如下：
[0089]
[0090] 式中： 和 分别表示第p个所述最优分解模式up的平均值和第q个所述最优分解模式uq的平均值；
[0091]
[0092] 式中：S表示所述最优分解模式u的平方包络谱， 和 分别表示第p个所述最优分解模式的平方包络谱Sp的平均值和第q个所述最优分解模式的平方包络谱Sq的平均值；
[0093] 步骤S12，对所述筛选分解模式组合中的两个最优分解模式进行计算，得到第一包络谱峭度值和第二包络谱峭度值，舍弃包络谱峭度值较小的分解模式，得到K个输出分解模式；
[0094] 所述包络谱峭度的计算公式如下：
[0095]
[0096] 步骤S13，对K进行判断，若K≤Q，则执行S14，若K＞Q，则执行S4；
[0097] 步骤S14，得到K个输出分解模式，基于K个输出分解模式构建包络谱，基于所述包络谱对待测轴承进行故障诊断。
[0098] 具体地，搭建如图1所示的试验台，该试验台包括电机、行星齿轮箱、平行轴齿轮箱、转速传感器、待测轴承、制动器等部件组成。电机经行星齿轮箱、平行轴齿轮箱两级增速后驱动轴承内圈运动，轴承外圈固定在轴承座上，加速度传感器固定在轴承座正上方，测量轴承的振动信号。滚动轴承即为待测轴承，内置外圈裂纹与内圈裂纹复合故障，采样频率为20k Hz，时长为1s，得到的双通道振动信号及其对应的包络谱图如图3所示，由图3可以观察到，信号中存在大量的噪声，且能量集中在100Hz以下的低频频段。只能显示出明显的内圈故障，无法找到任何关于外圈故障的有效信息
[0099] 基于图3所示的振动信号，采用本发明的方法得到的最优分解模式M1，M2和M3及其包络谱如图4所示，与图3相比，本发明方法揭示了被噪声完全淹没的外圈故障特征，且提取出了清晰的内圈故障频率的边频带信息，进一步增强了内圈故障特征。
[0100] 基于图3所示的振动信号，采用现有的FMD方法得到的最优分解模式及其包络谱如图5所示，可以发现，在最优分解结果中，FMD仍局限于突出的轴承内圈故障特征，而轴承外圈故障特征未能被成功提取出来，并且分解模式高度重复。因此，在这种情况下，FMD无法诊断出复合故障。
[0101] 基于图3所示的振动信号，采用现有的MVMD方法得到的最优分解模式及其包络谱如图6所示，由图可以观察到，MVMD的大部分分解模式都增强了内圈故障特征，仅有通道2的M6分解模式提取出了微弱的外圈故障特征，但其信噪比明显比本发明的方法得到的分解模式更低。并且，MVMD方法同时产生了大量冗余的重复模式与非故障模式。图1‑图6中，fI表示轴承内圈故障特征频率，fO表示轴承外圈故障特征频率，fr表示转频。
[0102] 本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0103] 以上所述的实施例仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。
