CN115901248B - 一种齿轮箱故障特征提取方法 - Google Patents

Abstract

本发明公开了一种齿轮箱故障特征提取方法，一种基于ISSA‑VMD‑MCKD的方法，其中ISSA算法为改进的麻雀搜索算法，VMD算法为变分模态分解算法，MCKD算法为最大相关峭度解卷积算法。所述齿轮箱故障特征提取方法包括顺序采用基于ISSA的VMD算法选取最优模态、基于ISSA的MCKD算法加强振动信号以及倒平方包络谱提取信号故障特征的方法。本发明通过对麻雀搜索算法的改进，以及基于ISSA对传统依靠经验确定的VMD和MCKD参数进行全局寻优，从而实现齿轮箱故障诊断正确率更高、效率更高、故障特征提取更充分的目的。

Description

一种齿轮箱故障特征提取方法

技术领域

本发明涉及齿轮箱故障特征提取技术领域，具体涉及一种齿轮箱故障特征提取方法。

背景技术

齿轮箱属于旋转机械类别。一般旋转机械类故障受原始数据噪声、采样信号非线性和不平稳状态等多类型交织因素的影响，故障特征提取不充分，导致诊断正确率和效率都不高。国内外研究采用小波降噪、经验模态分解(EMD)与改进的局部Fisher判别分析相结合、基于能量聚集度经验小波变换、基于小波包分析与BP神经网络、基于最大相关峭度解卷积(MCKD)和变分模态分解(VMD)、基于平方包络谱、基于麻雀搜索算法(SSA)优化支持向量机等故障诊断方法，虽然获得了一定的效果，但针对上述自适应选择问题还很难做到真正自适应，另一方面EMD、EEMD、LMD等属于递归模态分解，尚缺乏严格的数学理论。因此相关研究至今仍是一个热门问题。

齿轮箱广泛应用于工程机械等大型复杂机械装备中，低速重载以及恶劣工作环境经常导致齿轮箱关键部位出现严重故障，再加上强背景噪声的影响，大量有效振动信号被淹没在噪声中，给故障检测带来困难。另外，在实际齿轮箱的振动信号中，由于多个齿轮产生了多种转速和啮合频率，而且常常受到多个调制源的联合作用，形成了非对称的边带结构，检测信号功率谱中间包含很多大小和周期都不相同的频率结构，这就导致检测信号中混杂的周期分量很难被分辨出来。

当前，在齿轮箱故障特征提取技术领域，现有的传统故障诊断理论与技术不能有效解决齿轮箱所面临的诸多棘手问题，例如齿轮箱中振动传输路径复杂以及多模式混淆导致故障响应微弱、载荷大范围瞬时波动引起振动的强烈非平稳性、多对齿轮啮合的振动相互耦合造成振动明显的非线性，混杂的周期变量很难被分辨、噪声污染严重等问题。主要表现为故障特征提取不充分、诊断正确率和效率不高等现象。

随着物联网和信息技术的发展，对齿轮箱故障检测而言，传感器技术、检测技术、信号传输技术等都有了长足发展，信号处理计算能力也得到了显著进步，关键是缺少有效实用的故障信号分析处理的方法。

发明内容

有鉴于此，本发明提出了一种齿轮箱故障特征提取方法，通过基于ISSA(改进麻雀搜索算法)的VMD和MCKD算法解决齿轮箱中故障信号微弱、大量有效信号被背景噪声淹没以及故障信号所具有的线性和非线性特征等问题，将处理后的故障信号通过倒平方包络谱分析方法，更好地提取故障特征频率，识别相应的故障类型。，从而有效解决上述现有技术问题。

本发明设计的一种齿轮箱故障特征提取方法，所述方法基于ISSA-VMD-MCKD算法，也称为“一种基于ISSA-VMD-MCKD的齿轮箱故障特征提取方法”，其中ISSA算法为改进的麻雀搜索算法，VMD算法为变分模态分解算法，MCKD算法为最大相关峭度解卷积算法，其特征在于，所述故障特征提取方法包括如下步骤：

S1：采集齿轮箱运行时的振动信号；

S2：基于ISSA优化VMD算法，获取所述振动信号的最优模态；

S3：基于ISSA优化MCKD算法和所述最优模态，获取加强后的振动信号，特别是微弱振动信号获得明显加强；

S4：基于倒平方包络谱方法对所述加强后的振动信号提取振动信号故障特征；

所述ISSA的改进点包括在基本麻雀搜索算法的加入者位置更新中用前一时刻的最优位置替换传统正态分布随机数；所述基于ISSA优化VMD算法，包括对所述VMD算法信号中的分解层数K、惩罚因子a参数利用ISSA算法进行全局寻优的方法；所述基于ISSA优化MCKD算法，包括对所述MCKD算法中的滤波长度L、解卷积周期T参数利用ISSA算法进行全局寻优的方法。

其中，变分模态分解(VMD)是一种非递归式模态分解算法，它可以有效避免经验模态分解(EMD)等信号分解产生的模态混叠、端点效应等问题。而最大相关峭度解卷积(MCKD)通过解卷积运算突出被噪声淹没的连续冲击脉冲，提高原始信号的相关峭度值，该方法非常适用于提取微弱故障信号的连续瞬态冲击。采用两种方法组合对齿轮箱故障特征提取具有明显的优势，可以充分发挥VMD在降噪方面的优越性以及MCKD能够突出被噪声所掩盖的微弱脉冲的特性。该技术为成熟技术，可以参考文献“基于MCKD和VMD的滚动轴承微弱故障特征提取”(《振动与冲击》2017年，Vol.36，No.20)等。

本发明考虑到，利用VMD和MCKD虽然能够有效提取出齿轮箱中有效故障特征，但是这两种算法需要进行人为的设置一些参数，并且这些参数的设置对算法结果的影响非常大。具体而言，VMD算法需要设置惩罚因子a和模态分解层数K；MCKD需要提前设置滤波长度参数L和解卷积周期T。为了提高诊断效果，本发明提出了采用改进麻雀搜索算法(ISSA)分别对VMD和MCKD需要确定的参数进行全局寻优的方法，可减少人为不确定因素的影响，并有效提高算法的执行效率。

关于VMD和MCKD执行的先后顺序，传统方法一般是先利用MCKD进行微弱信号加强，再利用VMD进行模态分解。本发明考虑到若先采用MCKD算法可能会导致原始信号中的大量微弱噪声信号得到加强，有效信号会再次被噪声信号淹没，如此再进行VMD分解，VMD的分解效果会受到影响。所以本发明在设计过程中提出新方法之一就是首先考虑先进行VMD分解，然后对分解后的信号进行MCKD反卷积，保留最具特征的信号。另一方面同时也保留了传统的先采用MCKD算法，得到解卷积的信号后，再对该信号进行VMD分解的方案。从实验看先进行VMD模态分解后再进行MCKD反卷积微弱信号加强效果更好一些。

关于倒平方包络谱的使用：由于频谱图中故障特征频率信号一般较小，而包络谱则对冲击事件的故障比较敏感，包络谱图中故障特征频率的幅值较高，容易辨认。因此，相对于频谱分析，包络谱分析可剔除部分频率干扰。另外考虑到齿轮箱中每对齿轮啮合都将产生边频带，这就导致边频带会交叉分布在一起，仅靠包络谱分析还是很难清晰的观察特征频率，特别是能够凸显故障的特征频率。本发明提出了倒平方包络谱分析方法来提取信号的故障特征，倒平方包络谱能够较好地检测出频谱上的周期成分，并将原频谱上成族的边频带谱线简化成为单根谱线，便于观察，而齿轮发生故障时的振动频谱具有的边频带一般都具有等间隔的结构，利用倒平方包络谱这个优点，可以检测出功率谱中难以辨识的周期性信号。另外若只分析包络谱信号，故障特征频率聚集在低频区就很难分辨，不能直观的判断出齿轮箱的故障类型，而倒平方包络谱分析方法可以突出聚集在低频段的故障特征频率。经过倒平方包络谱分析方法后，能够更好地识别相应的故障类型。

所述VMD算法和MCKD算法均为现有技术，具体方法步骤可以参考文献“基于MCKD和VMD的滚动轴承微弱故障特征提取”(《振动与冲击》2017年，Vol.36，No.20)，其他类似文献均有相关描述。

所述VMD算法是一种完全非递归模式的信号分解方法，根据预设分解层数K，通过最优变分模型的迭代计算，将采集到的复杂数字信号x(t)分解为K个离散的模态u_k(t)，其中，1≤k≤K，t为时间，由于K可以预先设定，根据模态分解原理，如果K取值恰当，就可以有效的抑制模态混叠现象；

其中，A_k(t)为第k个分量的振幅，为第k个分量的相位；每个分量的中心频率为u_k(t)的瞬时频率；在不影响符号表达意义的情况下，为简便起见u_k(t)有时也记为u_k，ω_k(t)简记为ω_k；

所述VMD算法包括如下步骤：

S21：初始化λ¹＝0，n＝0，k＝0；

S22：n＝n+1(本表述为计算机语句)，开始整个算法的循环；

S23：k＝k+1(本表述为计算机语句)，直到k＝K，循环更新

S24：更新λ：

其中，Γ表示噪声容限，当信号含有强噪声时，可设定Γ＝0达到更好地去噪效果；

S25：判断是否满足收敛条件，如果满足收敛条件则停止迭代，否则返回步骤S22；所述收敛条件包括：

和/或循环次数限制；其中，ε为一个大于0的整数，代表精度；表示矩阵或行列式的2范数的平方，上标2为平方运算，下标2为2范数；所述收敛条件还可以包括迭代次数不超过某一上限，以防止迭代无限制循环。

其中，步骤S23所述循环更新包括如下分步骤：

S231：通过希尔伯特变换计算每个模态的解析信号，从而获得模态分量的单边频谱

其中：δ(t)为脉冲函数；t为时间；j为虚数单位；“*”表示卷积；

S232：通过单边频谱与算子进行频率混合，将各模态的中心带调制到相应的基带

S233：以分量相加等于原信号作为约束条件，构建约束变分模型如下：

约束条件为

其中，代表分解后的模态分量，代表各分量的中心频率；表示函数对t求偏导。

S234：为求解所述变分模型的最优解，引入拉格朗日乘子λ(t)及二次惩罚因子a，将约束性性变分问题变为非约束性变分问题；其中，λ(t)能够保证约束条件的严格性，a能够保证高斯噪声环境下信号重构的准确性；扩展的拉格朗日表达式为：

S235：利用交替方向乘子算法(ADDM)不断迭代求解如下自适应中心频率及各模态分量：

其中，分别是对应于x(t),λ(t)的傅里叶变换；还可以采用迭代自适应求解方法。

所述MCKD算法的基本原理为：

齿轮箱的振动信号一般为时间序列下的振动幅值，齿轮箱发生故障时因出现局部碰撞会产生周期性信号y_i。冲击信号y_i传递到传感器时不仅会因传输路径的影响而衰减，而且会掺杂大量的噪声，所以传感器采集到的信号可以表示为：

x_i＝hy_i+e

其中：x_i为传感器实际采集到的信号，h为原始信号在传输时经过的系统传递函数，e为齿轮箱工作环境下的噪声。

MCKD算法实质上是通过寻找一系列FIR滤波器f，使得原始周期性冲击序列的相关峭度最大，即寻找滤波器使实际信号x_i恢复到y_i，即

y＝f*x

其中：y，x分别为y_i，x_i的向量形式，向量长度为N，即i＝1,2,3,…,N；f＝[f₁,f₂,…,f_L]^T，当y的相关峭度达到最大时，则认为此时的f为期望的目标；此处*为卷积；L为滤波器长度参数；l∈[0,L]。

MCKD算法中，周期冲击信号y_i的相关峭度定义为：

其中：T为冲击信号周期；M为移位数；m∈[0,M]；增加M，可以增加算法序列脉冲，M过大会影响精度，一般M＝min{8,b/T}，其中b为采样点数。

对实际信号进行滤波，使得相关峭度最大，即

求解上式，等价于求解下列方程

可求得滤波器的最终系数，并表示为矩阵形式：

其中： r＝[0,T,…,mT]，上标T表示矩阵转置；

所述MCKD算法包括如下步骤：

S31：初始化解卷积周期T、位移数M及滤波器长度L等参数；

S32：根据MCKD算法原理计算输入信号x的X₀、X_T、

S33：计算滤波后的输出信号y；

S34：根据y计算β和Ψ；

S35：更新滤波器f的系数；

S36：若滤波前、后信号的相关峭度值ΔCK_M(T)小于阈值和/或迭代次数超过限制，结束迭代，否则重复步骤S33～S35。

进一步的，所述ISSA算法基本原理步骤包括：

麻雀搜索算法(SSA)的基本原理是：麻雀种群觅食时，其个体分为发现者和加入者两种角色；发现者承担搜索事物丰富的区域的责任，加入者能利用发现者获取食物丰富区域的信息；麻雀个体所对应的适应度值的好坏由能量储备的高低决定。

SA1：发现者位置更新包括如下式的方法：

式中：为第i个麻雀在第j维中的位置信息，麻雀总数为I；j＝1,2,…,d,其中d代表解的维度或待优化的参数个数；iter_max表示最大迭代次数；t表示现阶段的迭代次数；α表示均匀分布在(0,1]中的随机数；Q表示服从标准正态分布的随机数；T表示为1*d阶矩阵，且矩阵元素都为1；R₂表示预警值，取值大小在[0,1]范围内；ST表示安全值，取值处于[0,1]的范围内；当R₂<ST时，这意味着此时的觅食环境周围没有捕食者，发现者可以执行广泛的搜索操作；当R₂≥ST时，这表示种群中的一些麻雀已经发现了捕食者，并向种群中其他麻雀发出了警报，此时所有麻雀都需要迅速飞到其他安全的地方进行觅食；

SA2：加入者位置更新，如果发现危险，个体会向其他麻雀发出警告；当报警值大于安全值时，加入者会跟随发现者的路径到达安全区域；加入者位置更新包括如下式的方法：

其中：是前一时刻的全局最优位置，替换传统正态分布随机数；表示在d维空间中种群在进行第t次迭代后种群处于的最差位置；表示在第d维空间中种群第t+1次迭代后种群处于的最优位置；A为1*d的矩阵，其中每个元素随机赋值为1或-1，A⁺＝A^T(AA^T)^-1为A的伪逆矩阵，|·|表示取模运算；

当i>n/2时，传统公式为Q为服从正态分布的随机数，本发明将Q替换为前一时刻的最优位置根据前一时刻的最优位置，来寻找下一时刻的位置，可减少不必要的位置更新。

发现者和加入者时刻处于动态平衡之中，每只个体均在两种角色之间相互转换，但是两种角色个体数量占整个种群数量的比重是不变的；

加入者的适应度值越低，所在区域食物数量越少，位置越差；

加入者能够发现处于食物最丰富位置的发现者，从而获取该区域的食物或者围绕着发现者四周觅食；另一方面，发起者时时刻刻受到加入者的监控，对食物资源发起争夺。

SA3：危险来临时位置更新，危险来临时种群个体立即从种群边缘位置向安全区域迅速移动，位于种群中间的个体则会随机移动，向其他个体靠近，公式如下：

其中，β为服从标准正态分布的随机数，用于控制个体移动的步长；G表示个体移动的方向，G∈[-1,1]是一个随机数；ε是控制步长参数的一个极小常数，避免式中分母为零，其取值范围为[-0.5,0.5]；f_i表示第i只个体的适应度函数值，f_g和f_w分别用于表示麻雀种群的最优和最差适应度值f_i、f_g和f_w在求解过程中是变化的。

所述麻雀搜索算法为现有技术，具体可见参考文献“基于麻雀搜索算法优化支持向量机的滚动轴承故障诊断”(《科学技术与工程》2021年，Vol.21，No.10)，或文献“改进麻雀搜索算法在分布式电源配置中的应用”(《计算机工程与应用》2021年，Vol.57，No.20)，本发明的改进点与文献改进点不同，目的也不同。

进一步的，所述ISSA算法改进点还包括基于Levy飞行策略的加入者位置更新方法改进，包括用以下步骤SA21代替步骤SA2，所述步骤SA21包括：加入者位置更新如下式，

当i≤n/2时，现有技术基于Levy飞行策略的位置更新公式为： 本发明认为麻雀搜索算法出现陷入局部极值的问题，主要是因为算法在前期进行搜索时出现问题，因此对i≤n/2时的公式进行进一步优化；而i>n/2时依然用前一时刻的最优位置替换传统的Q。

式中：表示现阶段发现者的最优位置，α表示Levy飞行的步长，α>0，Levy(·)表示麻雀个体的随机搜寻轨迹，服从下式所示的Levy分布：

其中，β在(0,2)范围内的随机数，μ服从N(0,σ²)的正态分布，v服从N(0,1)标准正态分布，σ计算方法如下式：

其中：Γ是Gamma分布函数，用于模拟Levy分布；

麻雀搜索算法在搜索前期，麻雀个体位置变化较大，Levy飞行的步长取值应该较大，以达到快速定位最优解的目的；随着搜索进行到后期，麻雀个体位置变化较小，Levy飞行步长应该逐渐减小，以提高算法求解的精度。

进一步的，本发明在引入Levy飞行的基础上使用自适应动态步长优化麻雀搜索算法，自适应动态步长公式如下：

式中：c为现阶段的迭代次数；c_max为最大迭代次数。α表示步长大小，|·|表示取模运算，而传统的麻雀搜索算法通常设置α＝1。

麻雀搜索算法具有收敛速度快、误差小等优点，然而也会有陷入局部极值的问题。为了提高算法的全局搜索能力，降低出现局部最优解的可能性，本发明针对基于Levy飞行策略的麻雀搜索算法(也称为一种改进的麻雀搜索算法ISSA，有文献也简写为LSSA)，对SSA加入者位置更新公式中进行进一步优化，增加麻雀个体搜索时的扰动变异概率，使麻雀搜索算法在搜索后期具有较好种群多样性，提高其跳出局部最优的能力。基于Levy飞行策略的麻雀搜索算法为现有技术，具体可参考文献“改进麻雀搜索算法在分布式电源配置中的应用”(《计算机工程与应用》2021年，Vol.57，No.20)。

进一步的，定义无量纲指标包络谱峰值因子Ec为适应度函数，该指标同时考虑了信号冲击成分的周期性和强度，记信号包络谱幅值序列为X(z)(z＝1,2,3,…,Z，Z为最大序列数，则Ec表达为

所述基于ISSA优化VMD算法包括：将所述VMD算法信号的分解层数K、惩罚因子a两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入VMD算法获取输入信号的最优模态；

所述基于ISSA优化MCKD算法包括：将所述MCKD算法的滤波长度L、解卷积周期T两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入MCKD算法获取加强的振动信号；

所述步骤S2中将原始信号经过VMD算法分解后，基于包络谱峰值因子找出最优模态，然后将该最优模态作为步骤S3中MCKD算法的输入；

所述基于包络谱峰值因子判断是否满足优化要求和基于包络谱峰值因子找出最优模态的方法包括：判断包络谱峰值因子是否最小；

VMD算法中的惩罚因子a和模态分解层数K以及MCKD算法中的滤波长度参数L和解卷积周期T都需要进行人为的设置参数，且这些参数的取值对算法性能的影响非常大，由于此类故障诊断方法一般诊断一次都需要数分钟到数小时不等，若经验不足的人很有可能设置出错误参数，这就会极大降低算法寻优的性能，程序调试时间也会极大地延长。

这就需要我们通过其他算法来实现以上参数的自适应选择，麻雀搜索算法是近年来提出的一种智能优化算法，相比于类似于粒子群，遗传算法等优化方法，SSA算法具有可调参数少，易操作，鲁棒性强，应用范围广等优点。但是又考虑到SSA算法容易陷入到局部最优解，所以本发明在SSA的基础上进行修改，提高其跳出局部最优的能力，然后用改进的SSA算法开对VMD和MCKD算法进行参数优化。

进一步的，所述基于倒平方包络谱方法提取信号故障特征包括如下步骤：

S41：针对ISSA优化VMD-MCKD算法处理后的信号y(k)，其Hilbert变换为：

H[y(k)]＝F^-1{F{y(k)}h(k)}

其中：k＝1,2,3,…,L，L为信号采样数或信号长度，一般与采样频率以及采样时间相关；F{·}与F^-1{·}分别表示傅里叶变换与傅里叶逆变换，h(k)定义为

S42：构造解析信号

z(k)＝y(k)+jH[y(k)]

S43：平方包络信号y_SE(k)定义为解析信号z(k)模数的平方，即

y_SE(k)＝|z(k)|²

S44：平方包络谱Y_SE(n)则由平方包络信号y_SE(k)傅里叶变换幅值的平方得到，即

Y_SE(n)＝|F{y_SE(k)}|²

对平方包络信号Y_SE(n)作倒频谱分析即可得到倒平方包络谱。

S45：用C_p(n)来表示平方包络信号Y_SE(n)的功率谱函数的倒频谱，有

C_p(n)＝F^-1{log₁₀[Y_SE(n)]}

S46：倒平方包络谱为C_a(n)为

C_a(n)＝|F^-1{log₁₀[Y_SE(n)]}|

其中，F^-1{}为傅里叶逆变换，|·|表示取模运算。

一般频谱图中故障特征频率信号较小，而包络谱对冲击事件的故障比较敏感，包络谱图中故障特征频率的幅值较高，容易辨认。因此，相对于频谱分析，包络谱分析剔除了不必要的频率干扰。但考虑到齿轮箱中每对齿轮啮合都将产生边频带，这就导致边频带会交叉分布在一起，仅靠包络谱分析还是无法清晰的观察特征频率，特别是要凸显故障的特征频率。本发明提出了倒平方包络谱分析方法来提取信号的故障特征，倒平方包络谱能够较好地检测出频谱上的周期成分，将原来谱上成族的边频带谱线简化成为单根谱线，便于观察，而齿轮箱发生故障时的振动频谱具有的边频带一般都具有等间隔的结构，利用倒平方包络谱这个优点，可以检测出功率谱中难以辨识的周期性信号。即若只通过包络分析后得到的包络谱信号，故障特征频率聚集在低频区就很难分辨，不能直观的判断出轴承的故障类型，为了进一步对包络谱信号进行处理，本发明提出的倒平方包络谱分析方法可以突出聚集在低频段的故障特征频率。经过倒平方包络谱分析方法后，能够更好地识别相应的故障类型。

进一步的，所述故障类型包括但不局限于旋转机械齿轮箱齿轮缺齿、断齿、磨损、齿根裂纹、轴承内圈故障、外圈、滚珠、混合故障等。

另一方面，一种齿轮箱故障特征提取方法，所述方法基于ISSA-MCKD-VMD算法，也称为“一种基于ISSA-MCKD-VMD的齿轮箱故障特征提取方法”，其中ISSA算法为改进的麻雀搜索算法，VMD算法为变分模态分解算法，MCKD算法为最大相关峭度解卷积算法，其特征在于，所述故障特征提取方法包括如下步骤：

B1：采集齿轮箱运行时的振动信号；

B2：基于ISSA优化MCKD算法，获取加强后的振动信号；

B3：基于ISSA优化VMD算法和所述加强后的振动信号，获取所述振动信号的最优模态；

B4：基于倒平方包络谱方法对最优模态提取振动信号故障特征；

所述ISSA的改进点包括：在基本麻雀搜索算法的加入者位置更新中用前一时刻的最优位置替换传统正态分布随机数；所述基于ISSA优化VMD算法，包括对所述VMD算法信号中的分解层数K、惩罚因子a参数利用ISSA算法进行全局寻优的方法；所述基于ISSA优化MCKD算法，包括对所述MCKD算法中的滤波长度L、解卷积周期T参数利用ISSA算法进行全局寻优的方法。

其中，变分模态分解(VMD)是一种非递归式模态分解算法，它可以有效避免经验模态分解(EMD)等信号分解产生的模态混叠、端点效应等问题。而最大相关峭度解卷积(MCKD)通过解卷积运算突出被噪声淹没的连续冲击脉冲，提高原始信号的相关峭度值，该方法非常适用于提取微弱故障信号的连续瞬态冲击。采用两种方法组合对齿轮箱故障特征提取具有明显的优势，可以充分发挥VMD在降噪方面的优越性以及MCKD能够突出被噪声所掩盖的微弱脉冲的特性。该技术为成熟技术，可以参考文献“基于MCKD和VMD的滚动轴承微弱故障特征提取”(《振动与冲击》2017年，Vol.36，No.20)等。

本发明考虑到，利用VMD和MCKD虽然能够有效提取出齿轮箱中有效故障特征，但是这两种算法需要进行人为的设置一些参数，并且这些参数的设置对算法结果的影响非常大。具体而言，VMD算法需要设置惩罚因子a和模态分解层数K；MCKD需要提前设置滤波长度参数L和解卷积周期T。为了提高诊断效果，本发明提出了采用改进麻雀搜索算法(ISSA)分别对VMD和MCKD需要确定的参数进行全局寻优的方法，可减少人为不确定因素的影响，并有效提高算法的执行效率。

关于倒平方包络谱的使用：由于频谱图中故障特征频率信号一般较小，而包络谱则对冲击事件的故障比较敏感，包络谱图中故障特征频率的幅值较高，容易辨认。因此，相对于频谱分析，包络谱分析可剔除部分频率干扰。另外考虑到齿轮箱中每对齿轮啮合都将产生边频带，这就导致边频带会交叉分布在一起，仅靠包络谱分析还是很难清晰的观察特征频率，特别是能够凸显故障的特征频率。本发明提出了倒平方包络谱分析方法来提取信号的故障特征，倒平方包络谱能够较好地检测出频谱上的周期成分，并将原频谱上成族的边频带谱线简化成为单根谱线，便于观察，而齿轮发生故障时的振动频谱具有的边频带一般都具有等间隔的结构，利用倒平方包络谱这个优点，可以检测出功率谱中难以辨识的周期性信号。另外若只分析包络谱信号，故障特征频率聚集在低频区就很难分辨，不能直观的判断出齿轮箱的故障类型，而倒平方包络谱分析方法可以突出聚集在低频段的故障特征频率。经过倒平方包络谱分析方法后，能够更好地识别相应的故障类型。

所述VMD算法和MCKD算法均为现有技术，具体方法步骤可以参考文献“基于MCKD和VMD的滚动轴承微弱故障特征提取”(《振动与冲击》2017年，Vol.36，No.20)，其他类似文献均有相关描述。

所述VMD算法是一种完全非递归模式的信号分解方法，根据预设分解层数K，通过最优变分模型的迭代计算，将采集到的复杂数字信号x(t)分解为K个离散的模态u_k(t)，其中，1≤k≤K，t为时间，由于K可以预先设定，根据模态分解原理，如果K取值恰当，就可以有效的抑制模态混叠现象；

其中，A_k(t)为第k个分量的振幅，为第k个分量的相位；每个分量的中心频率为u_k(t)的瞬时频率；在不影响符号表达意义的情况下，为简便起见u_k(t)有时也记为u_k，ω_k(t)简记为ω_k；

所述VMD算法包括如下步骤：

B31：初始化λ¹＝0，n＝0，k＝0；

B32：n＝n+1(本表述为计算机语句)，开始整个算法的循环；

B33：k＝k+1(本表述为计算机语句)，直到k＝K，循环更新

B34：更新λ：

其中，Γ表示噪声容限，当信号含有强噪声时，可设定Γ＝0达到更好地去噪效果；

B35：判断是否满足收敛条件，如果满足收敛条件则停止迭代，否则返回步骤B32；所述收敛条件包括：

和/或循环次数限制；其中，ε为一个大于0的整数，代表精度；表示矩阵或行列式的2范数的平方，上标2为平方运算，下标2为2范数；所述收敛条件还可以包括迭代次数不超过某一上限，以防止迭代无限制循环。

其中，步骤B33所述循环更新包括如下分步骤：

B331：通过希尔伯特变换计算每个模态的解析信号，从而获得模态分量的单边频谱

其中：δ(t)为脉冲函数；t为时间；j为虚数单位；“*”表示卷积；

B332：通过单边频谱与算子进行频率混合，将各模态的中心带调制到相应的基带

B333：以分量相加等于原信号作为约束条件，构建约束变分模型如下：

约束条件为

其中，代表分解后的模态分量，代表各分量的中心频率；表示函数对t求偏导。

B334：为求解所述变分模型的最优解，引入拉格朗日乘子λ(t)及二次惩罚因子a，将约束性性变分问题变为非约束性变分问题；其中，λ(t)能够保证约束条件的严格性，a能够保证高斯噪声环境下信号重构的准确性；扩展的拉格朗日表达式为：

B335：利用交替方向乘子算法(ADDM)不断迭代求解如下自适应中心频率及各模态分量：

其中，分别是对应于x(t),λ(t)的傅里叶变换；还可以采用迭代自适应求解方法。

所述MCKD算法的基本原理为：

齿轮箱的振动信号一般为时间序列下的振动幅值，齿轮箱发生故障时因出现局部碰撞会产生周期性信号y_i。冲击信号y_i传递到传感器时不仅会因传输路径的影响而衰减，而且会掺杂大量的噪声，所以传感器采集到的信号可以表示为：

x_i＝hy_i+e

其中：x_i为传感器实际采集到的信号，h为原始信号在传输时经过的系统传递函数，e为齿轮箱工作环境下的噪声。

MCKD算法实质上是通过寻找一系列FIR滤波器f，使得原始周期性冲击序列的相关峭度最大，即寻找滤波器使实际信号x_i恢复到y_i，即

y＝f*x

其中：y，x分别为y_i，x_i的向量形式，向量长度为N，即i＝1,2,3,…,N；f＝[f₁,f₂,…,f_L]^T，当y的相关峭度达到最大时，则认为此时的f为期望的目标；此处*为卷积；L为滤波器长度参数；l∈[0,L]。

MCKD算法中，周期冲击信号y_i的相关峭度定义为：

其中：T为冲击信号周期；M为移位数；m∈[0,M]；增加M，可以增加算法序列脉冲，M过大会影响精度，一般M＝min{8,b/T}，其中b为采样点数。

对实际信号进行滤波，使得相关峭度最大，即

求解上式，等价于求解下列方程

可求得滤波器的最终系数，并表示为矩阵形式：

其中： r＝[0,T,…,mT]，上标T表示矩阵转置；

所述MCKD算法包括如下步骤：

B21：初始化解卷积周期T、位移数M及滤波器长度L等参数；

B22：根据MCKD算法原理计算输入信号x的X₀、X_T、

B23：计算滤波后的输出信号y；

B24：根据y计算β和Ψ；

B25：更新滤波器f的系数；

B26：若滤波前、后信号的相关峭度值ΔCK_M(T)小于阈值和/或迭代次数超过限制，结束迭代，否则重复步骤B23～B25。

进一步的，所述ISSA算法基本原理步骤包括：

麻雀搜索算法(SSA)的基本原理是：麻雀种群觅食时，其个体分为发现者和加入者两种角色；发现者承担搜索事物丰富的区域的责任，加入者能利用发现者获取食物丰富区域的信息；麻雀个体所对应的适应度值的好坏由能量储备的高低决定。

BA1：发现者位置更新包括如下式的方法：

式中：为第i个麻雀在第j维中的位置信息，麻雀总数为I；j＝1,2,…,d,其中d代表解的维度或待优化的参数个数；iter_max表示最大迭代次数；t表示现阶段的迭代次数；α表示均匀分布在(0,1]中的随机数；Q表示服从标准正态分布的随机数；T表示为1*d阶矩阵，且矩阵元素都为1；R₂表示预警值，取值大小在[0,1]范围内；ST表示安全值，取值处于[0,1]的范围内；当R₂<ST时，这意味着此时的觅食环境周围没有捕食者，发现者可以执行广泛的搜索操作；当R₂≥ST时，这表示种群中的一些麻雀已经发现了捕食者，并向种群中其他麻雀发出了警报，此时所有麻雀都需要迅速飞到其他安全的地方进行觅食；

BA2：加入者位置更新，如果发现危险，个体会向其他麻雀发出警告；当报警值大于安全值时，加入者会跟随发现者的路径到达安全区域；加入者位置更新包括如下式的方法：

其中：是前一时刻的全局最优位置，替换传统正态分布随机数；表示在d维空间中种群在进行第t次迭代后种群处于的最差位置；表示在第d维空间中种群第t+1次迭代后种群处于的最优位置；A为1*d的矩阵，其中每个元素随机赋值为1或-1，A⁺＝A^T(AA^T)^-1为A的伪逆矩阵，|·|表示取模运算；

当i>n/2时，传统公式为Q为服从正态分布的随机数，本发明将Q替换为前一时刻的最优位置根据前一时刻的最优位置，来寻找下一时刻的位置，可减少不必要的位置更新。

发现者和加入者时刻处于动态平衡之中，每只个体均在两种角色之间相互转换，但是两种角色个体数量占整个种群数量的比重是不变的；

加入者的适应度值越低，所在区域食物数量越少，位置越差；

加入者能够发现处于食物最丰富位置的发现者，从而获取该区域的食物或者围绕着发现者四周觅食；另一方面，发起者时时刻刻受到加入者的监控，对食物资源发起争夺。

BA3：危险来临时位置更新，危险来临时种群个体立即从种群边缘位置向安全区域迅速移动，位于种群中间的个体则会随机移动，向其他个体靠近，公式如下：

其中，β为服从标准正态分布的随机数，用于控制个体移动的步长；G表示个体移动的方向，G∈[-1,1]是一个随机数；ε是控制步长参数的一个极小常数，避免式中分母为零，其取值范围为[-0.5,0.5]；f_i表示第i只个体的适应度函数值，f_g和f_w分别用于表示麻雀种群的最优和最差适应度值f_i、f_g和f_w在求解过程中是变化的。

所述麻雀搜索算法为现有技术，具体可见参考文献“基于麻雀搜索算法优化支持向量机的滚动轴承故障诊断”(《科学技术与工程》2021年，Vol.21，No.10)，或文献“改进麻雀搜索算法在分布式电源配置中的应用”(《计算机工程与应用》2021年，Vol.57，No.20)，本发明的改进点与文献改进点不同，目的也不同。

进一步的，所述ISSA算法改进点还包括基于Levy飞行策略的加入者位置更新方法改进，包括用以下步骤BA21代替步骤BA2，所述步骤BA21包括：加入者位置更新如下式，

当i≤n/2时，现有技术基于Levy飞行策略的位置更新公式为： 本发明认为麻雀搜索算法出现陷入局部极值的问题，主要是因为算法在前期进行搜索时出现问题，因此对i≤n/2时的公式进行进一步优化；而i>n/2时依然用前一时刻的最优位置替换传统的Q。

式中：表示现阶段发现者的最优位置，α表示Levy飞行的步长，α>0，Levy(·)表示麻雀个体的随机搜寻轨迹，服从下式所示的Levy分布：

其中，β在(0,2)范围内的随机数，μ服从N(0,σ²)的正态分布，v服从N(0,1)标准正态分布，σ计算方法如下式：

其中：Γ是Gamma分布函数，用于模拟Levy分布；

麻雀搜索算法在搜索前期，麻雀个体位置变化较大，Levy飞行的步长取值应该较大，以达到快速定位最优解的目的；随着搜索进行到后期，麻雀个体位置变化较小，Levy飞行步长应该逐渐减小，以提高算法求解的精度。

进一步的，本发明在引入Levy飞行的基础上使用自适应动态步长优化麻雀搜索算法，自适应动态步长公式如下：

式中：c为现阶段的迭代次数；c_max为最大迭代次数。α表示步长大小，|·|表示取模运算，而传统的麻雀搜索算法通常设置α＝1。

麻雀搜索算法具有收敛速度快、误差小等优点，然而也会有陷入局部极值的问题。为了提高算法的全局搜索能力，降低出现局部最优解的可能性，本发明针对基于Levy飞行策略的麻雀搜索算法(也称为一种改进的麻雀搜索算法ISSA，有文献也简写为LSSA)，对SSA加入者位置更新公式中进行进一步优化，增加麻雀个体搜索时的扰动变异概率，使麻雀搜索算法在搜索后期具有较好种群多样性，提高其跳出局部最优的能力。基于Levy飞行策略的麻雀搜索算法为现有技术，具体可参考文献“改进麻雀搜索算法在分布式电源配置中的应用”(《计算机工程与应用》2021年，Vol.57，No.20)。

进一步的，定义无量纲指标包络谱峰值因子Ec为适应度函数，该指标同时考虑了信号冲击成分的周期性和强度，记信号包络谱幅值序列为X(z)(z＝1,2,3,…,Z，Z为最大序列数，则Ec表达为

所述基于ISSA优化VMD算法包括：将所述VMD算法信号的分解层数K、惩罚因子a两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入VMD算法获取输入信号的最优模态；

所述基于ISSA优化MCKD算法包括：将所述MCKD算法的滤波长度L、解卷积周期T两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入MCKD算法获取加强的振动信号；

所述步骤B3中将原始信号经过VMD算法分解后，基于包络谱峰值因子找出最优模态，然后将该最优模态作为步骤B2中MCKD算法的输入；

所述基于包络谱峰值因子判断是否满足优化要求和基于包络谱峰值因子找出最优模态的方法包括：判断包络谱峰值因子是否最小；

VMD算法中的惩罚因子a和模态分解层数K以及MCKD算法中的滤波长度参数L和解卷积周期T都需要进行人为的设置参数，且这些参数的取值对算法性能的影响非常大，由于此类故障诊断方法一般诊断一次都需要数分钟到数小时不等，若经验不足的人很有可能设置出错误参数，这就会极大降低算法寻优的性能，程序调试时间也会极大地延长。

这就需要我们通过其他算法来实现以上参数的自适应选择，麻雀搜索算法是近年来提出的一种智能优化算法，相比于类似于粒子群，遗传算法等优化方法，SSA算法具有可调参数少，易操作，鲁棒性强，应用范围广等优点。但是又考虑到SSA算法容易陷入到局部最优解，所以本发明在SSA的基础上进行修改，提高其跳出局部最优的能力，然后用改进的SSA算法开对VMD和MCKD算法进行参数优化。

进一步的，所述基于倒平方包络谱方法提取信号故障特征包括如下步骤：

B41：针对ISSA优化VMD-MCKD算法处理后的信号y(k)，其Hilbert变换为：

H[y(k)]＝F^-1{F{y(k)}h(k)}

其中：k＝1,2,3,…,L，L为信号采样数或信号长度，一般与采样频率以及采样时间相关；F{·}与F^-1{·}分别表示傅里叶变换与傅里叶逆变换，h(k)定义为

B42：构造解析信号

z(k)＝y(k)+jH[y(k)]

B43：平方包络信号y_SE(k)定义为解析信号z(k)模数的平方，即

y_SE(k)＝|z(k)|²

B44：平方包络谱Y_SE(n)则由平方包络信号y_SE(k)傅里叶变换幅值的平方得到，即

Y_SE(n)＝|F{y_SE(k)}|²

对平方包络信号Y_SE(n)作倒频谱分析即可得到倒平方包络谱。

B45：用C_p(n)来表示平方包络信号Y_SE(n)的功率谱函数的倒频谱，有

C_p(n)＝F^-1{log₁₀[Y_SE(n)]}

B46：倒平方包络谱为C_a(n)为

C_a(n)＝|F^-1{log₁₀[Y_SE(n)]}|

其中，F^-1{}为傅里叶逆变换，|·|表示取模运算。

一般频谱图中故障特征频率信号较小，而包络谱对冲击事件的故障比较敏感，包络谱图中故障特征频率的幅值较高，容易辨认。因此，相对于频谱分析，包络谱分析剔除了不必要的频率干扰。但考虑到齿轮箱中每对齿轮啮合都将产生边频带，这就导致边频带会交叉分布在一起，仅靠包络谱分析还是无法清晰的观察特征频率，特别是要凸显故障的特征频率。本发明提出了倒平方包络谱分析方法来提取信号的故障特征，倒平方包络谱能够较好地检测出频谱上的周期成分，将原来谱上成族的边频带谱线简化成为单根谱线，便于观察，而齿轮箱发生故障时的振动频谱具有的边频带一般都具有等间隔的结构，利用倒平方包络谱这个优点，可以检测出功率谱中难以辨识的周期性信号。即若只通过包络分析后得到的包络谱信号，故障特征频率聚集在低频区就很难分辨，不能直观的判断出轴承的故障类型，为了进一步对包络谱信号进行处理，本发明提出的倒平方包络谱分析方法可以突出聚集在低频段的故障特征频率。经过倒平方包络谱分析方法后，能够更好地识别相应的故障类型。

进一步的，所述故障类型包括但不局限于旋转机械齿轮箱齿轮缺齿、断齿、磨损、齿根裂纹、轴承内圈故障、外圈、滚珠、混合故障等。

另一方面，一种齿轮箱故障诊断系统，其特征在于：该系统执行如上所述的“一种基于ISSA-VMD-MCKD的齿轮箱故障特征提取方法”，或“一种基于ISSA-MCKD-VMD的齿轮箱故障特征提取方法”，实现齿轮箱故障特征提取或故障诊断的功能。

本发明的有益效果在于：通过对麻雀搜索算法的改进，一方面提高其运算效率，另一方面降低其陷入局部极值的概率；利用改进后的麻雀搜索算法，分别对VMD算法中的惩罚因子a和模态分解层数K以及MCKD算法中的滤波长度参数L和解卷积周期T参数进行全局寻优，改变了传统方法需要进行人为设置、不确定性较大的问题，可以更好地实现故障特征提取；针对齿轮故障时的振动频谱边频带一般都具有等间隔的结构特征，利用倒平方包络谱能够较好地检测频谱周期成分的优点，可以检测出功率谱中难以辨识的周期性信号，能够较好地克服传统包络谱信号分析或平方包络谱信号分析中，故障特征频率聚集在低频区难以分辨的问题。本发明提出的基于ISSA-VMD-MCKD和ISSA-MCKD-VMD的齿轮箱故障特征提取方法实现了齿轮箱故障诊断正确率更高、效率更高、故障特征提取更充分的目的。

附图说明

图1是一种基于ISSA-VMD-MCKD的齿轮箱故障特征提取方法原理框图；

图2是一种基于ISSA-MCKD-VMD的齿轮箱故障特征提取方法原理框图；

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

实施例1：一种基于ISSA-VMD-MCKD的齿轮箱故障特征提取方法

如图1所示，一种齿轮箱故障特征提取方法，所述方法基于ISSA-VMD-MCKD算法，其中ISSA算法为改进的麻雀搜索算法，VMD算法为变分模态分解算法，MCKD算法为最大相关峭度解卷积算法，其特征在于，所述故障特征提取方法包括如下步骤：

S1：采集齿轮箱运行时的振动信号；

S2：基于ISSA优化VMD算法，获取所述振动信号的最优模态；

S3：基于ISSA优化MCKD算法和所述最优模态，获取加强后的振动信号，特别是微弱振动信号获得明显加强；

S4：基于倒平方包络谱方法对所述加强后的振动信号提取振动信号故障特征；

所述ISSA的改进点包括在基本麻雀搜索算法的加入者位置更新中用前一时刻的最优位置替换传统正态分布随机数；所述基于ISSA优化VMD算法，包括对所述VMD算法信号中的分解层数K、惩罚因子a参数利用ISSA算法进行全局寻优的方法；所述基于ISSA优化MCKD算法，包括对所述MCKD算法中的滤波长度L、解卷积周期T参数利用ISSA算法进行全局寻优的方法。

其中，变分模态分解(VMD)是一种非递归式模态分解算法，它可以有效避免经验模态分解(EMD)等信号分解产生的模态混叠、端点效应等问题。而最大相关峭度解卷积(MCKD)通过解卷积运算突出被噪声淹没的连续冲击脉冲，提高原始信号的相关峭度值，该方法非常适用于提取微弱故障信号的连续瞬态冲击。采用两种方法组合对齿轮箱故障特征提取具有明显的优势，可以充分发挥VMD在降噪方面的优越性以及MCKD能够突出被噪声所掩盖的微弱脉冲的特性。该技术为成熟技术，可以参考文献“基于MCKD和VMD的滚动轴承微弱故障特征提取”(《振动与冲击》2017年，Vol.36，No.20)等。

所述VMD算法是一种完全非递归模式的信号分解方法，根据预设分解层数K，通过最优变分模型的迭代计算，将采集到的复杂数字信号x(t)分解为K个离散的模态u_k(t)，其中，1≤k≤K，t为时间，由于K可以预先设定，根据模态分解原理，如果K取值恰当，就可以有效的抑制模态混叠现象；

其中，A_k(t)为第k个分量的振幅，为第k个分量的相位；每个分量的中心频率为u_k(t)的瞬时频率；在不影响符号表达意义的情况下，为简便起见u_k(t)有时也记为u_k，ω_k(t)简记为ω_k；

所述VMD算法包括如下步骤：

S21：初始化λ¹＝0，n＝0，k＝0；本实施例根据被测齿轮箱经验模态进行初始化，所述经验模态一般可通过测量数据的简要分析即可得到，也可以根据转速、齿轮箱各齿轮的齿数及其啮合方式经理论分析得到，对于某一固定的齿轮箱每次故障分析初始化模态可以是固定的；

S22：n＝n+1(本表述为计算机语句)，开始整个算法的循环；

S23：k＝k+1(本表述为计算机语句)，直到k＝K，循环更新

S24：更新λ：

其中，Γ表示噪声容限，当信号含有强噪声时，可设定Γ＝0达到更好地去噪效果；

S25：判断是否满足收敛条件，如果满足收敛条件则停止迭代，否则返回步骤S22；所述收敛条件包括：

和/或循环次数限制；其中，ε为一个大于0的整数，代表精度；表示矩阵或行列式的2范数的平方，上标2为平方运算，下标2为2范数；所述收敛条件还可以包括迭代次数不超过某一上限，以防止迭代无限制循环。

其中，步骤S23所述循环更新包括如下分步骤：

S231：通过希尔伯特变换计算每个模态的解析信号，从而获得模态分量的单边频谱

其中：δ(t)为脉冲函数；t为时间；j为虚数单位；“*”表示卷积；

S232：通过单边频谱与算子进行频率混合，将各模态的中心带调制到相应的基带

S233：以分量相加等于原信号作为约束条件，构建约束变分模型如下：

约束条件为

其中，代表分解后的模态分量，代表各分量的中心频率；表示函数对t求偏导。

S234：为求解所述变分模型的最优解，引入拉格朗日乘子λ(t)及二次惩罚因子a，将约束性性变分问题变为非约束性变分问题；其中，λ(t)能够保证约束条件的严格性，a能够保证高斯噪声环境下信号重构的准确性；扩展的拉格朗日表达式为：

S235：利用交替方向乘子算法(ADDM)不断迭代求解如下自适应中心频率及各模态分量：

其中，分别是对应于x(t),λ(t)的傅里叶变换；还可以采用迭代自适应求解方法。

所述MCKD算法的基本原理为：

齿轮箱的振动信号一般为时间序列下的振动幅值，齿轮箱发生故障时因出现局部碰撞会产生周期性信号y_i。冲击信号y_i传递到传感器时不仅会因传输路径的影响而衰减，而且会掺杂大量的噪声，所以传感器采集到的信号可以表示为：

x_i＝hy_i+E

其中：x_i为传感器实际采集到的信号，h为原始信号在传输时经过的系统传递函数，e为齿轮箱工作环境下的噪声。

MCKD算法实质上是通过寻找一系列FIR滤波器f，使得原始周期性冲击序列的相关峭度最大，即寻找滤波器使实际信号x_i恢复到y_i，即

y＝f*x

其中：y，x分别为y_i，x_i的向量形式，向量长度为N，即i＝1,2,3,…,N；f＝[f₁,f₂,…,f_L]^T，当y的相关峭度达到最大时，则认为此时的f为期望的目标；此处*为卷积；L为滤波器长度参数；l∈[0,L]。

MCKD算法中，周期冲击信号y_i的相关峭度定义为：

其中：T为冲击信号周期；M为移位数；m∈[0,M]；增加M，可以增加算法序列脉冲，M过大会影响精度，一般M＝min{8,b/T}，其中b为采样点数。

对实际信号进行滤波，使得相关峭度最大，即

求解上式，等价于求解下列方程

可求得滤波器的最终系数，并表示为矩阵形式：

其中： r＝[0,T,…,mT]，上标T表示矩阵转置；

所述MCKD算法包括如下步骤：

S31：初始化解卷积周期T、位移数M及滤波器长度L等参数；

S32：根据MCKD算法原理计算输入信号x的X₀、X_T、

S33：计算滤波后的输出信号y；

S34：根据y计算β和Ψ；

S35：更新滤波器f的系数；

S36：若滤波前、后信号的相关峭度值ΔCK_M(T)小于阈值和/或迭代次数超过限制，结束迭代，否则重复步骤S33～S35。

优选的，所述ISSA算法基本原理步骤包括：

SA1：发现者位置更新包括如下式的方法：

式中：为第i个麻雀在第j维中的位置信息，麻雀总数为I；j＝1,2,…,d,其中d代表解的维度或待优化的参数个数，本实施例对VMD和MCKD各2个参数进行优化，因此取d＝2；iter_max表示最大迭代次数；t表示现阶段的迭代次数；α表示均匀分布在(0,1]中的随机数；Q表示服从标准正态分布的随机数；T表示为1*d阶矩阵，且矩阵元素都为1；R₂表示预警值，取值大小在[0,1]范围内，本实施例取计算机系统随机函数rand()生成的随机数；ST表示安全值，取值处于[0,1]的范围内，本实施例取ST＝0.6；当R₂<ST时，这意味着此时的觅食环境周围没有捕食者，发现者可以执行广泛的搜索操作；当R₂≥ST时，这表示种群中的一些麻雀已经发现了捕食者，并向种群中其他麻雀发出了警报，此时所有麻雀都需要迅速飞到其他安全的地方进行觅食；

SA2：加入者位置更新，如果发现危险，个体会向其他麻雀发出警告；当报警值大于安全值时，加入者会跟随发现者的路径到达安全区域；加入者位置更新包括如下式的方法：

其中：是前一时刻的全局最优位置，替换传统正态分布随机数；表示在d维空间中种群在进行第t次迭代后种群处于的最差位置；表示在第d维空间中种群第t+1次迭代后种群处于的最优位置；A为1*d的矩阵，其中每个元素随机赋值为1或-1，A⁺＝A^T(AA^T)^-1为A的伪逆矩阵，|·|表示取模运算；

SA3：危险来临时位置更新，危险来临时种群个体立即从种群边缘位置向安全区域迅速移动，位于种群中间的个体则会随机移动，向其他个体靠近，公式如下：

其中，β为服从标准正态分布的随机数，用于控制个体移动的步长；G表示个体移动的方向，G∈[-1,1]是一个随机数；ε是控制步长参数的一个极小常数，避免式中分母为零，其取值范围为[-0.5,0.5]；f_i表示第i只个体的适应度函数值，f_g和f_w分别用于表示麻雀种群的最优和最差适应度值f_i、f_g和f_w在求解过程中是变化的。

优选的，定义无量纲指标包络谱峰值因子Ec为适应度函数，该指标同时考虑了信号冲击成分的周期性和强度，记信号包络谱幅值序列为X(z)(z＝1,2,3,…,Z，Z为最大序列数，则Ec表达为

所述基于ISSA优化VMD算法包括：将所述VMD算法信号的分解层数K、惩罚因子a两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入VMD算法获取输入信号的最优模态；

所述基于ISSA优化MCKD算法包括：将所述MCKD算法的滤波长度L、解卷积周期T两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入MCKD算法获取加强的振动信号；

所述步骤S2中将原始信号经过VMD算法分解后，基于包络谱峰值因子找出最优模态，然后将该最优模态作为步骤S3中MCKD算法的输入；

所述基于包络谱峰值因子判断是否满足优化要求和基于包络谱峰值因子找出最优模态的方法包括：判断包络谱峰值因子是否最小，选取包络谱峰值因子最小值对应的参数组合为最优参数组合，或者对应的模态为最优模态。

优选的，所述基于倒平方包络谱方法提取信号故障特征包括如下步骤：

S41：针对ISSA优化VMD-MCKD算法处理后的信号y(k)，其Hilbert变换为：

H[y(k)]＝F^-1{F{y(k)}h(k)}

其中：k＝1,2,3,…,L，L为信号采样数或信号长度，一般与采样频率以及采样时间相关，本实施例为信号采样点数；F{·}与F^-1{·}分别表示傅里叶变换与傅里叶逆变换，h(k)定义为

S42：构造解析信号

z(k)＝y(k)+jH[y(k)]

S43：平方包络信号y_SE(k)定义为解析信号z(k)模数的平方，即

y_SE(k)＝|z(k)|²

S44：平方包络谱Y_SE(n)则由平方包络信号y_SE(k)傅里叶变换幅值的平方得到，即

Y_SE(n)＝|F{y_SE(k)}|²

对平方包络信号Y_SE(n)作倒频谱分析即可得到倒平方包络谱。

S45：用C_p(n)来表示平方包络信号Y_SE(n)的功率谱函数的倒频谱，有

C_p(n)＝F^-1{log₁₀[Y_SE(n)]}

S46：倒平方包络谱为C_a(n)为

C_a(n)＝|F^-1{log₁₀[Y_SE(n)]}|

其中，F^-1{}为傅里叶逆变换，|·|表示取模运算。

以下是本实施例“一种基于ISSA-VMD-MCKD的齿轮箱故障特征提取方法”的一个完整步骤描述：

S1：采集齿轮箱运行时的振动信号；

S2：基于ISSA优化VMD算法，获取所述振动信号的最优模态；包括：

S21：VMD参数初始化λ¹＝0，n＝0，k＝0，关于K、a取值，一般K的取值范围是[3,10]，a的取值范围是[1000,1500]，本实施例用计算机系统随机函数rand()进行[K,a]的初始化，即用公式[K,a]＝int(rand()*([max_K,max_a]-[min_K,min_a])+[min_K,min_a])计算得到；其中max_K、max_a、min_K、min_a分别为K、a的取值范围上限和下限，int()表示取整；

S22：n＝n+1(本表述为计算机语句)，开始整个算法的循环；

S23：k＝k+1(本表述为计算机语句)，直到k＝K，循环更新包括：

S231：通过希尔伯特变换计算每个模态的解析信号，从而获得模态分量的单边频谱；

S232：通过单边频谱与算子进行频率混合，将各模态的中心带调制到相应的基带；

S233：以分量相加等于原信号作为约束条件，构建约束变分模型；

S234：为求解所述变分模型的最优解，引入拉格朗日乘子；

S235：利用交替方向乘子算法(ADDM)不断迭代求解如下自适应中心频率及各模态分量；

S24：更新λ；

S25：判断是否满足收敛条件，不满足，回到S22循环；满足，循环结束；本实施例VMD算法内循环次数限制设为不超过100次；

SA0：根据VMD获取的模态分量求解包络谱后再计算适应度函数，本实施例适应度函数为包络谱峰值因子Ec；设置ISSA待优化量X＝[K,a]，初始化参数，包括解的维度d＝2，最大迭代次数iter_max＝100，麻雀总数量I＝100，等，执行ISSA，包括：

SA1：发现者位置更新；

SA2：加入者位置更新；

SA3：危险来临时位置更新；

SA4：获取待解量X＝[K,a]的优化更新值，根据迭代条件判断是否结束循环，所述迭代条件包括迭代次数限制和/或计算结果门限限制；不满足迭代结束条件，则用[K,a]的优化更新值回到步骤S22重新执行循环；满足迭代结束条件，对所有[K,a]组合值对应的最优模态的适应度函数，选取适应度函数最小的组合，获得其对应VMD算法最优模态作为步骤S3中MCKD算法的输入；所述最优模态的选取方法包括计算所有对应模态的适应度函数，选取适应度函数最小的模态作为最优模态的方法；

S3：基于ISSA优化MCKD算法和所述最优模态，获取加强后的振动信号，特别是微弱振动信号获得明显加强；包括：

S31：参数初始化，本实施例取位移数M＝3，解卷积周期T和滤波器长度L使用公式int((Valuemax-Valuemin)*rand()+Valuemin)进行计算，其中计算T时Valuemax＝150，Valuemin＝90；计算L时Valuemax＝1000，Valuemin＝10；int()为取整函数，rand()为计算机系统自带随机函数，可产生0～1的随机数；

S32：根据MCKD算法原理计算输入信号x的X₀、X_T、

S33：计算滤波后的输出信号y；

S34：根据y计算β和Ψ；

S35：更新滤波器f的系数；

S36：若滤波前、后信号的相关峭度值ΔCK_M(T)小于阈值和/或迭代次数超过限制，结束迭代，否则重复步骤S33～S35；本实施例迭代次数上限设为100次，相关峭度值变化阈值设为0.001；最小迭代次数10次；

SA0：根据MCKD获取的加强后的振动信号求解包络谱后再计算适应度函数，本实施例适应度函数为包络谱峰值因子Ec；设置ISSA待优化量X＝[T,L]，初始化参数，包括解的维度d＝2，最大迭代次数iter_max＝100，麻雀总数量I＝100，等，执行ISSA，包括：

SA1：发现者位置更新；

SA2：加入者位置更新；

SA3：危险来临时位置更新；

SA4：获取待解量X＝[T,L]的优化更新值，根据迭代条件判断是否结束循环，所述迭代条件包括迭代次数限制和/或计算结果门限限制；不满足迭代结束条件，则用[T,L]的优化更新值回到步骤S32重新执行循环；满足迭代结束条件，对所有[T,L]组合值对应的加强后的振动信号的适应度函数，选取适应度函数最小的组合，获得其对应MCKD算法最优加强后的振动信号作为步骤S4中基于倒平方包络谱方法的输入；所述最优加强后的振动信号的选取方法包括计算所有对应[T,L]组合的适应度函数，选取适应度函数最小的振动信号作为最优加强后的振动信号的方法；

S4：基于倒平方包络谱方法对所述加强后的振动信号提取振动信号故障特征，包括步骤S41至S46。

优选的，所述故障类型包括但不局限于旋转机械齿轮箱齿轮缺齿、断齿、磨损、齿根裂纹、轴承内圈故障、外圈、滚珠、混合故障等。

实施例2：另一种基于ISSA-VMD-MCKD的齿轮箱故障特征提取方法

与实施例1的区别在于ISSA的改进点不同。

优选的，所述ISSA算法改进点还包括基于Levy飞行策略的加入者位置更新方法改进，包括用以下步骤SA21代替步骤SA2，所述步骤SA21包括：加入者位置更新如下式，

式中：表示现阶段发现者的最优位置，α表示Levy飞行的步长，α>0，Levy(·)表示麻雀个体的随机搜寻轨迹，服从下式所示的Levy分布：

其中，β在(0,2)范围内的随机数，μ服从N(0,σ²)的正态分布，v服从N(0,1)标准正态分布，σ计算方法如下式：

其中：Γ是Gamma分布函数，用于模拟Levy分布；

优选的，本发明在引入Levy飞行的基础上使用自适应动态步长优化麻雀搜索算法，自适应动态步长公式如下：

式中：c为现阶段的迭代次数；c_max为最大迭代次数。α表示步长大小，|·|表示取模运算，而传统的麻雀搜索算法通常设置α＝1。

实施例3：一种基于ISSA-MCKD-VMD的齿轮箱故障特征提取方法

与实施例1的主要区别在于“基于ISSA优化MCKD算法”和“基于ISSA优化VMD算法”执行先后顺序调换。

如图2所示，另一种齿轮箱故障特征提取方法，所述方法基于ISSA-MCKD-VMD算法，也称为“一种基于ISSA-MCKD-VMD的齿轮箱故障特征提取方法”，其中ISSA算法为改进的麻雀搜索算法，VMD算法为变分模态分解算法，MCKD算法为最大相关峭度解卷积算法，其特征在于，所述故障特征提取方法包括如下步骤：

B1：采集齿轮箱运行时的振动信号；

B2：基于ISSA优化MCKD算法，获取加强后的振动信号；

B3：基于ISSA优化VMD算法和所述加强后的振动信号，获取所述振动信号的最优模态；

B4：基于倒平方包络谱方法对最优模态提取振动信号故障特征；

所述ISSA的改进点包括：在基本麻雀搜索算法的加入者位置更新中用前一时刻的最优位置替换传统正态分布随机数；所述基于ISSA优化VMD算法，包括对步骤S3所述的VMD算法信号中的分解层数K、惩罚因子a参数利用ISSA算法进行全局寻优的方法；所述基于ISSA优化MCKD算法，包括对步骤S2所述的MCKD算法中的滤波长度L、解卷积周期T参数利用ISSA算法进行全局寻优的方法。

其中，变分模态分解(VMD)是一种非递归式模态分解算法，它可以有效避免经验模态分解(EMD)等信号分解产生的模态混叠、端点效应等问题。而最大相关峭度解卷积(MCKD)通过解卷积运算突出被噪声淹没的连续冲击脉冲，提高原始信号的相关峭度值，该方法非常适用于提取微弱故障信号的连续瞬态冲击。采用两种方法组合对齿轮箱故障特征提取具有明显的优势，可以充分发挥VMD在降噪方面的优越性以及MCKD能够突出被噪声所掩盖的微弱脉冲的特性。该技术为成熟技术，可以参考文献“基于MCKD和VMD的滚动轴承微弱故障特征提取”(《振动与冲击》2017年，Vol.36，No.20)等。

所述VMD算法是一种完全非递归模式的信号分解方法，根据预设分解层数K，通过最优变分模型的迭代计算，将采集到的复杂数字信号x(t)分解为K个离散的模态u_k(t)，其中，1≤k≤K，t为时间，由于K可以预先设定，根据模态分解原理，如果K取值恰当，就可以有效的抑制模态混叠现象；

其中，A_k(t)为第k个分量的振幅，为第k个分量的相位；每个分量的中心频率为u_k(t)的瞬时频率；在不影响符号表达意义的情况下，为简便起见u_k(t)有时也记为u_k，ω_k(t)简记为ω_k；

所述VMD算法包括如下步骤：

B31：初始化λ¹＝0，n＝0，k＝0；

B32：n＝n+1(本表述为计算机语句)，开始整个算法的循环；

B33：k＝k+1(本表述为计算机语句)，直到k＝K，循环更新

B34：更新λ：

其中，Γ表示噪声容限，当信号含有强噪声时，可设定Γ＝0达到更好地去噪效果；

B35：判断是否满足收敛条件，如果满足收敛条件则停止迭代，否则返回步骤B32；所述收敛条件包括：

和/或循环次数限制；其中，ε为一个大于0的整数，代表精度；表示矩阵或行列式的2范数的平方，上标2为平方运算，下标2为2范数；所述收敛条件还可以包括迭代次数不超过某一上限，以防止迭代无限制循环。

其中，步骤B33所述循环更新包括如下分步骤：

B331：通过希尔伯特变换计算每个模态的解析信号，从而获得模态分量的单边频谱

其中：δ(t)为脉冲函数；t为时间；j为虚数单位；“*”表示卷积；

B332：通过单边频谱与算子进行频率混合，将各模态的中心带调制到相应的基带

B333：以分量相加等于原信号作为约束条件，构建约束变分模型如下：

约束条件为

其中，代表分解后的模态分量，代表各分量的中心频率；表示函数对t求偏导。

B334：为求解所述变分模型的最优解，引入拉格朗日乘子λ(t)及二次惩罚因子a，将约束性性变分问题变为非约束性变分问题；其中，λ(t)能够保证约束条件的严格性，a能够保证高斯噪声环境下信号重构的准确性；扩展的拉格朗日表达式为：

B335：利用交替方向乘子算法(ADDM)不断迭代求解如下自适应中心频率及各模态分量：

其中，分别是对应于x(t),λ(t)的傅里叶变换；还可以采用迭代自适应求解方法。

所述MCKD算法的基本原理为：

齿轮箱的振动信号一般为时间序列下的振动幅值，齿轮箱发生故障时因出现局部碰撞会产生周期性信号y_i。冲击信号y_i传递到传感器时不仅会因传输路径的影响而衰减，而且会掺杂大量的噪声，所以传感器采集到的信号可以表示为：

x_i＝hy_i+e

其中：x_i为传感器实际采集到的信号，h为原始信号在传输时经过的系统传递函数，e为齿轮箱工作环境下的噪声。

MCKD算法实质上是通过寻找一系列FIR滤波器f，使得原始周期性冲击序列的相关峭度最大，即寻找滤波器使实际信号x_i恢复到y_i，即

y＝f*x

其中：y，x分别为y_i，x_i的向量形式，向量长度为N，即i＝1,2,3,…,N；f＝[f₁,f₂,…,f_L]^T，当y的相关峭度达到最大时，则认为此时的f为期望的目标；此处*为卷积；L为滤波器长度参数；l∈[0,L]。

MCKD算法中，周期冲击信号y_i的相关峭度定义为：

其中：T为冲击信号周期；M为移位数；m∈[0,M]；增加M，可以增加算法序列脉冲，M过大会影响精度，一般M＝min{8,b/T}，其中b为采样点数。

对实际信号进行滤波，使得相关峭度最大，即

求解上式，等价于求解下列方程

可求得滤波器的最终系数，并表示为矩阵形式：

其中： r＝[0,T,…,mT]，上标T表示矩阵转置；

所述MCKD算法包括如下步骤：

B21：初始化解卷积周期T、位移数M及滤波器长度L等参数；

B22：根据MCKD算法原理计算输入信号x的X₀、X_T、

B23：计算滤波后的输出信号y；

B24：根据y计算β和Ψ；

B25：更新滤波器f的系数；

B26：若滤波前、后信号的相关峭度值ΔCK_M(T)小于阈值和/或迭代次数超过限制，结束迭代，否则重复步骤B23～B25。

优选的，所述ISSA算法基本原理步骤包括：

BA1：发现者位置更新包括如下式的方法：

式中：为第i个麻雀在第j维中的位置信息，麻雀总数为I；j＝1,2,…,d,其中d代表解的维度或待优化的参数个数；iter_max表示最大迭代次数；t表示现阶段的迭代次数；α表示均匀分布在(0,1]中的随机数；Q表示服从标准正态分布的随机数；T表示为1*d阶矩阵，且矩阵元素都为1；R₂表示预警值，取值大小在[0,1]范围内，本实施例取计算机系统随机函数rand()生成的随机数；ST表示安全值，取值处于[0,1]的范围内，本实施例取ST＝0.6；当R₂<ST时，这意味着此时的觅食环境周围没有捕食者，发现者可以执行广泛的搜索操作；当R₂≥ST时，这表示种群中的一些麻雀已经发现了捕食者，并向种群中其他麻雀发出了警报，此时所有麻雀都需要迅速飞到其他安全的地方进行觅食；

BA2：加入者位置更新，如果发现危险，个体会向其他麻雀发出警告；当报警值大于安全值时，加入者会跟随发现者的路径到达安全区域；加入者位置更新包括如下式的方法：

其中：是前一时刻的全局最优位置，替换传统正态分布随机数；表示在d维空间中种群在进行第t次迭代后种群处于的最差位置；表示在第d维空间中种群第t+1次迭代后种群处于的最优位置；A为1*d的矩阵，其中每个元素随机赋值为1或-1，A⁺＝A^T(AA^T)^-1为A的伪逆矩阵，|·|表示取模运算；

BA3：危险来临时位置更新，危险来临时种群个体立即从种群边缘位置向安全区域迅速移动，位于种群中间的个体则会随机移动，向其他个体靠近，公式如下：

其中，β为服从标准正态分布的随机数，用于控制个体移动的步长；G表示个体移动的方向，G∈[-1,1]是一个随机数；ε是控制步长参数的一个极小常数，避免式中分母为零，其取值范围为[-0.5,0.5]；f_i表示第i只个体的适应度函数值，f_g和f_w分别用于表示麻雀种群的最优和最差适应度值f_i、f_g和f_w在求解过程中是变化的。

优选的，定义无量纲指标包络谱峰值因子Ec为适应度函数，该指标同时考虑了信号冲击成分的周期性和强度，记信号包络谱幅值序列为X(z)(z＝1,2,3,…,Z，Z为最大序列数，则Ec表达为

所述基于ISSA优化VMD算法包括：将所述VMD算法信号的分解层数K、惩罚因子a两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入VMD算法获取输入信号的最优模态；

所述基于ISSA优化MCKD算法包括：将所述MCKD算法的滤波长度L、解卷积周期T两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入MCKD算法获取加强的振动信号；

所述步骤B3中将原始信号经过VMD算法分解后，基于包络谱峰值因子找出最优模态，然后将该最优模态作为步骤B2中MCKD算法的输入；

所述基于包络谱峰值因子判断是否满足优化要求和基于包络谱峰值因子找出最优模态的方法包括：判断包络谱峰值因子是否最小，选取包络谱峰值因子最小值对应的参数组合为最优参数组合，或者对应的模态为最优模态。

优选的，所述基于倒平方包络谱方法提取信号故障特征包括如下步骤：

B41：针对ISSA优化VMD-MCKD算法处理后的信号y(k)，其Hilbert变换为：

H[y(k)]＝F^-1{F{y(k)}h(k)}

其中：k＝1,2,3,…,L，L为信号采样数或信号长度，一般与采样频率以及采样时间相关；F{·}与F^-1{·}分别表示傅里叶变换与傅里叶逆变换，h(k)定义为

B42：构造解析信号

z(k)＝y(k)+jH[y(k)]

B43：平方包络信号y_SE(k)定义为解析信号z(k)模数的平方，即

y_SE(k)＝|z(k)|²

B44：平方包络谱Y_SE(n)则由平方包络信号y_SE(k)傅里叶变换幅值的平方得到，即

Y_SE(n)＝|F{y_SE(k)}|²

对平方包络信号Y_SE(n)作倒频谱分析即可得到倒平方包络谱。

B45：用C_p(n)来表示平方包络信号Y_SE(n)的功率谱函数的倒频谱，有

C_p(n)＝F^-1{log₁₀[Y_SE(n)]}

B46：倒平方包络谱为C_a(n)为

C_a(n)＝|F^-1{log₁₀[Y_SE(n)]}|

其中，F^-1{}为傅里叶逆变换，|·|表示取模运算。

以下是本实施例“一种基于ISSA-MCKD-VMD的齿轮箱故障特征提取方法”的一个完整步骤描述：

B1：采集齿轮箱运行时的振动信号；

B2：基于ISSA优化MCKD算法，获取加强后的振动信号；包括：

B21：参数初始化，本实施例取位移数M＝3，解卷积周期T和滤波器长度L使用公式int((Valuemax-Valuemin)*rand()+Valuemin)进行计算，其中计算T时Valuemax＝150，Valuemin＝90；计算L时Valuemax＝1000，Valuemin＝10；int()为取整函数，rand()为计算机系统自带随机函数，可产生0～1的随机数；

B22：根据MCKD算法原理计算输入信号x的X₀、X_T、

B23：计算滤波后的输出信号y；

B24：根据y计算β和Ψ；

B25：更新滤波器f的系数；

B26：若滤波前、后信号的相关峭度值ΔCK_M(T)小于阈值和/或迭代次数超过限制，结束迭代，否则重复步骤B23～B25；本实施例迭代次数上限设为100次，相关峭度值变化阈值设为0.001；最小迭代次数10次；

BA0：根据MCKD获取的加强后的振动信号求解包络谱后再计算适应度函数，本实施例适应度函数为包络谱峰值因子Ec；设置ISSA待优化量X＝[T,L]，初始化参数，包括解的维度d＝2，最大迭代次数iter_max＝100，麻雀总数量I＝100，等，执行ISSA，包括：

BA1：发现者位置更新；

BA2：加入者位置更新；

BA3：危险来临时位置更新；

BA4：获取待解量X＝[T,L]的优化更新值，根据迭代条件判断是否结束循环，所述迭代条件包括迭代次数限制和/或计算结果门限限制；不满足迭代结束条件，则用[T,L]的优化更新值回到步骤B22重新执行循环；满足迭代结束条件，对所有[T,L]组合值对应的加强后的振动信号的适应度函数，选取适应度函数最小的组合，获得其对应MCKD算法最优加强后的振动信号作为步骤B3中VMD算法的输入；所述最优加强后的振动信号的选取方法包括计算所有对应[T,L]组合的适应度函数，选取适应度函数最小的振动信号作为最优加强后的振动信号的方法；

B3：基于ISSA优化VMD算法和所述加强后的振动信号，获取所述振动信号的最优模态；

B31：VMD参数初始化λ¹＝0，n＝0，k＝0，关于K、a取值，一般K的取值范围是[3,10]，a的取值范围是[1000,1500]，本实施例用计算机系统随机函数rand()进行[K,a]的初始化，即用公式[K,a]＝int(rand()*([max_K,max_a]-[min_K,min_a])+[min_K,min_a])计算得到；其中max_K、max_a、min_K、min_a分别为K、a的取值范围上限和下限，int()表示取整；

B32：n＝n+1(本表述为计算机语句)，开始整个算法的循环；

B33：k＝k+1(本表述为计算机语句)，直到k＝K，循环更新包括：

B331：通过希尔伯特变换计算每个模态的解析信号，从而获得模态分量的单边频谱；

B332：通过单边频谱与算子进行频率混合，将各模态的中心带调制到相应的基带；

B333：以分量相加等于原信号作为约束条件，构建约束变分模型；

B334：为求解所述变分模型的最优解，引入拉格朗日乘子；

B335：利用交替方向乘子算法(ADDM)不断迭代求解如下自适应中心频率及各模态分量；

B34：更新λ；

B35：判断是否满足收敛条件，不满足，回到B32循环；满足，循环结束；本实施例VMD算法内循环次数限制设为不超过100次；

BA0：根据VMD获取的模态分量求解包络谱后再计算适应度函数，本实施例适应度函数为包络谱峰值因子Ec；设置ISSA待优化量X＝[K,a]，初始化参数，包括解的维度d＝2，最大迭代次数iter_max＝100，麻雀总数量I＝100，等，执行ISSA，包括：

BA1：发现者位置更新；

BA2：加入者位置更新；

BA3：危险来临时位置更新；

BA4：获取待解量X＝[K,a]的优化更新值，根据迭代条件判断是否结束循环，所述迭代条件包括迭代次数限制和/或计算结果门限限制；不满足迭代结束条件，则用[K,a]的优化更新值回到步骤B32重新执行循环；满足迭代结束条件，对所有[K,a]组合值对应的最优模态的适应度函数，选取适应度函数最小的组合，获得其对应VMD算法最优模态作为步骤B4中基于倒平方包络谱方法的输入；所述最优模态的选取方法包括计算所有对应模态的适应度函数，选取适应度函数最小的模态作为最优模态的方法；

B4：基于倒平方包络谱方法对所述加强后的振动信号提取振动信号故障特征，包括步骤B41至B46。

优选的，所述故障类型包括但不局限于旋转机械齿轮箱齿轮缺齿、断齿、磨损、齿根裂纹、轴承内圈故障、外圈、滚珠、混合故障等。

实施例4：另一种基于ISSA-MCKD-VMD的齿轮箱故障特征提取方法

与实施例3的区别在于ISSA的改进点不同。

优选的，所述ISSA算法改进点还包括基于Levy飞行策略的加入者位置更新方法改进，包括用以下步骤BA21代替步骤BA2，所述步骤BA21包括：加入者位置更新如下式，

式中：表示现阶段发现者的最优位置，α表示Levy飞行的步长，α>0，Levy(·)表示麻雀个体的随机搜寻轨迹，服从下式所示的Levy分布：

其中，β在(0,2)范围内的随机数，μ服从N(0,σ²)的正态分布，v服从N(0,1)标准正态分布，σ计算方法如下式：

其中：Γ是Gamma分布函数，用于模拟Levy分布；

优选的，本发明在引入Levy飞行的基础上使用自适应动态步长优化麻雀搜索算法，自适应动态步长公式如下：

式中：c为现阶段的迭代次数；c_max为最大迭代次数。α表示步长大小，|·|表示取模运算，而传统的麻雀搜索算法通常设置α＝1。

实施例5：一种齿轮箱故障特征提取系统

一种齿轮箱故障特征提取系统，该系统执行实施例1、2、3、4中任意一个或多个所述的方法。

本发明的基本原理在发明内容中已经在对应的主要改进点部分做了相应阐述，这里不再重述。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和润饰，包括优选方法的组合，如对麻雀搜索算法改进点重新进行取舍、对VMD和MCKD算法进行取舍等，这些改进和润饰也应视为本发明的保护范围。

Claims (10)

1.一种齿轮箱故障特征提取方法，所述方法基于ISSA-VMD-MCKD算法，其中ISSA算法为改进的麻雀搜索算法，VMD算法为变分模态分解算法，MCKD算法为最大相关峭度解卷积算法，其特征在于，所述故障特征提取方法包括如下步骤：

S1：采集齿轮箱运行时的振动信号；

S2：基于ISSA优化VMD算法，获取所述振动信号的最优模态；

S3：基于ISSA优化MCKD算法和所述最优模态，获取加强后的振动信号；

S4：基于倒平方包络谱方法对所述加强后的振动信号提取振动信号故障特征；

所述ISSA的改进点包括在基本麻雀搜索算法的加入者位置更新中用前一时刻的最优位置替换传统正态分布随机数；所述基于ISSA优化VMD算法，包括对所述VMD算法信号中的分解层数K、惩罚因子a参数利用ISSA算法进行全局寻优的方法；所述基于ISSA优化MCKD算法，包括对所述MCKD算法中的滤波长度L、解卷积周期T参数利用ISSA算法进行全局寻优的方法。

2.根据权利要求1所述的一种齿轮箱故障特征提取方法，其特征在于，所述ISSA算法步骤包括：

SA1：发现者位置更新包括如下式的方法：

式中：为第i个麻雀在第j维中的位置信息，麻雀总数为I；j＝1,2,…,d,其中d代表解的维度或待优化的参数个数；iter_max表示最大迭代次数；t表示现阶段的迭代次数；α表示均匀分布在(0,1]中的随机数；Q表示服从标准正态分布的随机数；T表示为1*d阶矩阵，且矩阵元素都为1；R₂表示预警值，取值大小在[0,1]范围内；ST表示安全值，取值处于[0,1]的范围内；当R₂<ST时，这意味着此时的觅食环境周围没有捕食者，发现者可以执行广泛的搜索操作；当R₂≥ST时，这表示种群中的一些麻雀已经发现了捕食者，并向种群中其他麻雀发出了警报，此时所有麻雀都需要迅速飞到其他安全的地方进行觅食；

SA2：加入者位置更新，如果发现危险，个体会向其他麻雀发出警告；当报警值大于安全值时，加入者会跟随发现者的路径到达安全区域；加入者位置更新包括如下式的方法：

其中：是前一时刻的全局最优位置，替换传统正态分布随机数；表示在d维空间中种群在进行第t次迭代后种群处于的最差位置；表示在第d维空间中种群第t+1次迭代后种群处于的最优位置；A为1*d的矩阵，其中每个元素随机赋值为1或-1，A⁺＝A^T(AA^T)^-1为A的伪逆矩阵，|·|表示取模运算；

SA3：危险来临时位置更新，危险来临时种群个体立即从种群边缘位置向安全区域迅速移动，位于种群中间的个体则会随机移动，向其他个体靠近，公式如下：

其中，β为服从标准正态分布的随机数，用于控制个体移动的步长；G表示个体移动的方向，G∈[-1,1]是一个随机数；ε是控制步长参数的一个极小常数，避免式中分母为零，其取值范围为[-0.5,0.5]；f_i表示第i只个体的适应度函数值，f_g和f_w分别用于表示麻雀种群的最优和最差适应度值f_i、f_g和f_w在求解过程中是变化的。

3.根据权利要求2所述的一种齿轮箱故障特征提取方法，其特征在于，所述ISSA算法改进点还包括基于Levy飞行策略的加入者位置更新方法改进，包括用以下步骤SA21代替步骤SA2，所述步骤SA21包括：加入者位置更新如下式，

式中：表示现阶段发现者的最优位置，α表示Levy飞行的步长，α>0，Levy(·)表示麻雀个体的随机搜寻轨迹，服从下式所示的Levy分布：

其中，β在(0,2)范围内的随机数，μ服从N(0,σ²)的正态分布，v服从N(0,1)标准正态分布，σ计算方法如下式：

其中：Γ是Gamma分布函数，用于模拟Levy分布；

式中：c为现阶段的迭代次数；c_max为最大迭代次数,α表示步长大小，|·|表示取模运算。

4.根据权利要求1、2、3中任意一项所述的一种齿轮箱故障特征提取方法，其特征在于，定义无量纲指标包络谱峰值因子Ec为适应度函数，记信号包络谱幅值序列为X(z)(z＝1,2,3,…,Z)，Z为最大序列数，则Ec表达为

所述基于ISSA优化VMD算法包括：将所述VMD算法信号的分解层数K、惩罚因子a两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入VMD算法获取输入信号的最优模态；

所述基于ISSA优化MCKD算法包括：将所述MCKD算法的滤波长度L、解卷积周期T两个参数组合成向量作为ISSA算法的待优化位置量，基于包络谱峰值因子判断是否满足优化要求，获取优化值后再代入MCKD算法获取加强的振动信号；

所述步骤S2中将原始信号经过VMD算法分解后，基于包络谱峰值因子找出最优模态，然后将该最优模态作为步骤S3中MCKD算法的输入；

所述基于包络谱峰值因子判断是否满足优化要求和基于包络谱峰值因子找出最优模态的方法包括：判断包络谱峰值因子是否最小。

5.根据权利要求4所述的一种齿轮箱故障特征提取方法，其特征在于，所述基于倒平方包络谱方法提取信号故障特征包括如下步骤：

S41：针对ISSA优化VMD-MCKD算法处理后的信号y(k)，其Hilbert变换为：

H[y(k)]＝F^-1{F{y(k)}h(k)}

其中：k＝1,2,3,…,L，L为信号采样数或信号长度；F{·}与F^-1{·}分别表示傅里叶变换与傅里叶逆变换，h(k)定义为

S42：构造解析信号

z(k)＝y(k)+jH[y(k)]

S43：平方包络信号y_SE(k)定义为解析信号z(k)模数的平方，即

y_SE(k)＝|z(k)|²

S44：平方包络谱Y_SE(n)则由平方包络信号y_SE(k)傅里叶变换幅值的平方得到，即

Y_SE(n)＝|F{y_SE(k)}|²

S45：用C_p(n)来表示平方包络信号Y_SE(n)的功率谱函数的倒频谱，有

C_p(n)＝F^-1{log₁₀[Y_SE(n)]}

S46：倒平方包络谱为C_a(n)为

C_a(n)＝|F^-1{log₁₀[Y_SE(n)]}|

其中，F^-1{}为傅里叶逆变换，|·|表示取模运算。

6.根据权利要求1、2、3、5中任意一项所述的一种齿轮箱故障特征提取方法，其特征在于，故障类型包括旋转机械齿轮箱齿轮缺齿、断齿、磨损、齿根裂纹、轴承内圈故障、外圈、滚珠、混合故障。

7.根据权利要求4所述的一种齿轮箱故障特征提取方法，其特征在于，故障类型包括旋转机械齿轮箱齿轮缺齿、断齿、磨损、齿根裂纹、轴承内圈故障、外圈、滚珠、混合故障。

8.一种齿轮箱故障特征提取方法，所述方法基于ISSA-MCKD-VMD算法，其中ISSA算法为改进的麻雀搜索算法，VMD算法为变分模态分解算法，MCKD算法为最大相关峭度解卷积算法，其特征在于，所述故障特征提取方法包括如下步骤：

B1：采集齿轮箱运行时的振动信号；

B2：基于ISSA优化MCKD算法，获取加强后的振动信号；

B3：基于ISSA优化VMD算法和所述加强后的振动信号，获取所述振动信号的最优模态；

B4：基于倒平方包络谱方法对最优模态提取振动信号故障特征；

所述ISSA的改进点包括：在基本麻雀搜索算法的加入者位置更新中用前一时刻的最优位置替换传统正态分布随机数；所述基于ISSA优化VMD算法，包括对所述VMD算法信号中的分解层数K、惩罚因子a参数利用ISSA算法进行全局寻优的方法；所述基于ISSA优化MCKD算法，包括对所述MCKD算法中的滤波长度L、解卷积周期T参数利用ISSA算法进行全局寻优的方法。

9.根据权利要求8所述的一种齿轮箱故障特征提取方法，其特征在于，故障类型包括旋转机械齿轮箱齿轮缺齿、断齿、磨损、齿根裂纹、轴承内圈故障、外圈、滚珠、混合故障。

10.一种齿轮箱故障特征提取系统，其特征在于：该系统执行如权利要求1～9中任一权利要求所述的方法。