CN111597651A - 一种基于hwpso-svdd模型的滚动轴承性能退化评估方法 - Google Patents

Abstract

本发明涉及一种基于HWPSO‑SVDD模型的滚动轴承性能退化评估方法，包括以下步骤：S1：获取正常滚动轴承振动信号的特征指标，通过PCA算法得到融合特征指标，建立训练样本；S2：利用HWPSO优化算法对SVDD模型的惩罚参数C和核参数δ进行寻优；S3：利用寻优得到的惩罚参数C和核参数δ，构建HWPSO‑SVDD模型；S4：将训练样本输入HWPSO‑SVDD模型进行模型训练，得到训练完成的性能退化评估模型；S5：获取待测滚动轴承振动信号的融合特征指标，利用训练完成的性能退化评估模型，计算得到待测滚动轴承的性能退化指标V_c值；S6：设置自适应报警阈值，当V_c值超过自适应报警阈值时，对轴承的早期退化状态做出预警，本发明具有提高预测准确性和鲁棒性等优点。

Description

一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法

技术领域

本发明涉及滚动轴承性能评估领域，尤其是涉及一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法。

背景技术

机械设备应用于人类生活、工作和生产的方方面面，并在其中扮演了举足轻重的角色。目前，机械设备正朝着大型化、高速化、精密化、系统化、连续化和自动化方向发展，机械设备的运行环境越来越复杂多变，更是为设备的健康管理提出了新的挑战。随着设备的运行，零件老化、可靠性降低、剩余寿命减少等问题渐渐显现，设备能否持续安全高效地工作、维护能否及时有效地执行，亟待人类解决。设备一旦失效，将会直接影响整个机械系统的正常工作，甚至引起发恶性事件，造成生命和财产的严重损失。因此，设备维修体制正由传统的定期维修或事后维修向基于状态的视情维修转变，而作为建立合理维修策略的前提，设备性能退化评估也开始备受关注。

滚动轴承作为旋转机械中的关键零部件之一，其性能状态的好坏直接影响整台设备的运行可靠性。一般来说，滚动轴承在使用过程中都会经历从正常到退化直至失效的过程，而这期间通常要经历一系列不同的性能退化状态。如果能够在这个过程中对轴承的退化程度进行定量评估，则可以使设备维护策略的制定具有针对性。为其设计与维护提供指导依据，避免人员伤亡和设备的损失。

目前，工程上常用时域指标或频域指标来监测滚动轴承的运行状态。传统时域和频域特征中的一些特征指标，例如均方根、峭度等以及频域的傅里叶变换方法等虽然也能作为特征指标，但是单个的特征指标包含的信息量过于单一，不能全面的反映滚动轴承内部蕴含的丰富信息，所以必须要构建一个全面且有效的特征指标来反映滚动轴承的退化趋势。

支持向量数据描述(SVDD)是在支持向量机(SVM)基础上发展起来的一种有效的单值分类方法，旨在寻找包容目标类的最优超球体，使得被描述对象尽可能多地包容在这个超球体内。相比于其他的评估模型，SVDD只需正常样本进行模型训练。然而目前SVDD模型的参数优化容易陷入局部最优的缺陷，或者由于人为选择模型参数的盲目性，而导致模型的预测结果不准确。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，能够有效提高滚动轴承退化趋势评估的准确性。

本发明的目的可以通过以下技术方案来实现：

一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，包括以下步骤：

S1：获取正常滚动轴承振动信号的特征指标，通过PCA算法得到融合特征指标，建立训练样本；

S2：利用HWPSO优化算法对SVDD模型的惩罚参数C和核参数δ进行寻优；

S3：利用寻优得到的惩罚参数C和核参数δ，构建HWPSO-SVDD模型；

S4：将训练样本输入HWPSO-SVDD模型进行模型训练，得到训练完成的性能退化评估模型；

S5：获取待测滚动轴承振动信号的融合特征指标，利用训练完成的性能退化评估模型，计算得到待测滚动轴承的性能退化指标V_c值；

S6：设置自适应报警阈值线，根据V_c值对轴承的早期退化状态做出及时准确的预警。

进一步地，所述的滚动轴承振动信号的特征指标包括时域特征指标和频域特征指标，所述的时域特征指标包括有量纲时域特征指标和无量纲时域特征指标，所述的频域特征指标包括中心频率、重心频率、频率方差、频域幅值平均值和均方根率。

进一步地，所述的步骤S1具体包括：

S11：构建特征指标矩阵Y：

其中，y_ij为第j个特征指标中的第i个数据，n为每一种特征指标中的包含的数据个数，D为特征指标的种类数；

S12：计算特征指标矩阵Y中各特征指标的均值：

其中，

为第j个特征指标中数据均值；

S13：求解特征指标的方差

S14：对特征指标矩阵Y标准化，得到标准特征指标矩阵Y₀：

S15：计算特征指标协方差矩阵S：

S16：获取特征指标协方差矩阵S对角线上的各特征值λ₁，λ₂，...，λ_n，并从大到小依次排列；

S17：计算各主成分的贡献率：

其中，其中为w_r第r个主成分对应的贡献率，λ_i为第i个特征值，L为所构造的主成分数。

S18：选择贡献率高于80％的主成分，作为融合特征指标。

进一步地，所述的步骤S2具体包括：

S21：建立适应度函数F_sv；

S22：初始化设置改进粒子群算法参数和鲸鱼群算法参数，并分别设置PSOgbest适应度和领头鲸适应度为无穷大；

S23：通过改进粒子群算法和鲸鱼群算法的迭代，对SVDD模型的惩罚参数C和核参数δ寻优，得到最优解。

进一步优选地，所述的适应度函数F_sv的表达式为：

其中，N_sv为SVDD模型里最小化支持向量个数，G为训练样本数。

更进一步地，所述的步骤S23具体包括：

S231：计算惯性权重w、认知加速因子c₁、社会加速度因子c₂：

其中，w_min为惯性权重的最小值，其值为0.4，w_max为惯性权重的最大值，其值为0.9，iter_max为HWPSO最大迭代次数，iter为当前HWPSO迭代次数，c_1i为认知加速度因子的初始值，c_1f为认知加速度因子的最终值，c_2i为社会加速度因子的初始值，c_2f为社会加速度因子的最终值，c_1i和c_2f的值为2.5，c_1f和c_2i的值为0.5；

S232：获取每个粒子的当前适应度，判断粒子当前适应度是否小于该粒子pbest适应度，若是，则用该粒子当前适应度及其位置替换该粒子pbest适应度及其位置；

S233：判断当前最优粒子pbest适应度是否小于PSOgbest适应度，若是，则用该粒子pbest适应度替换PSOgbest适应度；

S234：更新每个粒子的速度和位置；

S235：根据边界检查修正更新全局最佳位置，若修正的最新全局最佳位置越界，则用最大值或最小值代替；

S236：利用粒子种群的最佳位置初始化鲸鱼种群算法，并通过鲸鱼种群算法的迭代，对粒子种群最佳位置PSO gbest进行优化；

S237：判断当前HWPSO迭代次数iter是否达到HWPSO最大迭代次数iter_max，若是，则执行步骤S238，否则当前HWPSO迭代次数iter的值加1，并返回执行步骤S231；

S238：输出粒子种群最佳位置PSOgbest，其对应的坐标即是SVDD模型最优的惩罚参数C和核参数δ。

更进一步地，其特征在于，所述的步骤S236具体包括：

S2361：判断gbest适应度是否小于领头鲸适应度，若是，则将gbest适应度及其位置设定为领头鲸适应度和位置；

S2362：计算WOA最大迭代次数im：

im＝[A×(iter)+B]

其中，A、B为常数；

S2363：获取每条鲸鱼的当前适应度，判断鲸鱼当前适应度是否小于领头鲸适应度，若是，则将该鲸鱼当前适应度及其位置设定为领头鲸的适应度和位置；

S2364：更新当前鲸鱼群个体的空间位置；

S2365：根据边界检查修正更新领头鲸位置，若修正的领头鲸位置越界，则用最大值或最小值代替；

S2366：判断当前WOA迭代次数t是否达到WOA最大迭代次数im，若是，则执行步骤S237，否则当前WOA迭代次数t的值加1，并返回执行步骤S2363。

更进一步地，所述的更新每个粒子的速度和位置的计算式为：

其中，

为粒子i在第k次迭代中第d维的当前速度，

为粒子i在第k次迭代中第d维的当前位置，

为粒子i在第k次迭代中第d维的更新后速度，即粒子i在第k+1次迭代中第d维的当前速度，

为粒子i在第k次迭代中第d维的更新后位置，即粒子i在第k+1次迭代中第d维的当前位置，rand₁为0到1之间随机数，pbest_ia为粒子i第d维的局部最佳位置，rand₂为0到1之间随机数，gbest为全局最佳位置；

所述的更新当前鲸鱼群个体的空间位置，具体包括：

若概率因子p_m＜0.5，且系数向量A₂满足|A₂|≥1时，则更新鲸鱼群个体空间位置的表达式为：

X(t+1)＝X_rand(t)-A₂D₁

D₁＝|C₂X_rand(t)-X(t)|

若概率因子p_m＜0.5，且系数向量A₂满足|A₂|＜1时，则更新鲸鱼群个体空间位置的表达式为：

X(t+1)＝X^*(t)-A₂D₂

D₂＝|C₂X^*(t)-X(t)|

C₂＝2r

A₂＝2ar-a

其中，概率因子p_m为0到1的随机数，X(t+1)为更新的鲸鱼个体位置，X(t)为当前鲸群个体空间位置，X_rand(t)为当前鲸群中鲸鱼个体随机位置，X^*(t)为迄今鲸鱼个体最佳位置，t为当前WOA迭代次数，A₂和C₂为系数向量，D₁表示个体X(t)与随机位置X_rand(t)的距离D₂表示包围步长，a为常数，并在迭代过程中从2减少到0，a_min＝0，a_max＝2，r为取值范围为[0，1]的随机向量；

若概率因子p_m≥0.5，则更新鲸鱼群个体空间位置的表达式为：

X(t+1)＝D₃e^blcos(2πl)+X^*(t)

D₃＝|X^*(t)-X(t)|

l＝(a₂-1)×rand+1

其中，b为常数，其值为1，l为介于[-1，1]之间的随机数，a₂为常数，D₃表示鲸群第i头鲸目前最佳位置到猎物之间的距离，rand为0到1间的随机数。

进一步地，所述的HWPSO-SVDD模型的目标函数为：

s.t.(z_i-c)^T(z_i-c)≤R²+ε_i，ε_i≥0，i＝1，2，…，M

其中，z_i为目标样本，M为训练样本个数，C为惩罚因子，ε_i为松弛变量，c为超球面中心，R为超球面半径，K(·)为核函数，α_i为拉格朗日乘子，z_sv为支持向量(位于超球面上的样本称为支持向量)。

更进一步地，所述的待测滚动轴承的性能退化指标V_c值的计算式为：

其中，V_c为融合特征指标距离超球体中心的距离，z_q为待测样本

V_c是一个连续变化的参数，表示设备偏离正常状态的程度。设定报警阈值，能够对轴承的健康状态进行实时的监控。根据统计学中的3σ法则，对于一个均值为

方差为σ²的高斯随机变量而言，取值落在区间

内的概率为99.73％。一旦某个取值超过该范围，有理由认为该取值来自其他变量。同样可以假设处于相近性能退化状态下的V_c值也符合一定的正态分布，一旦连续多个V_c值超过前面的V_c值界定的3σ取值范围(单个V_c值的超限可能是受外界干扰影响)，则认为轴承的性能退化状态已发生较大改变。该阈值不需要事先设定，根据实际所测得的V_c值动态改变，是自适应的。由于V_c是性能退化程度的增函数，所以只考虑阈值上限，即

自适应报警阈值的计算公式为：

其中，t_s表示可以确定的轴承处于正常状态的时刻，T(t)表示t时刻的V_c值，mean、std分别表示求均值和标准偏差。阈值T_h的求解分为3个阶段：第1阶段数据来源于早期无故障状态，是一个固定值；第2阶段将t时刻的T(t)与t-1时刻的T_h(t-1)进行比较，若T(t)在T_h(t-1)范围内，则将T(t)纳入原始数据，计算T_h(t)；第3阶段为其后连续N_u个V_c值均超限，则定义t＝t_e为轴承性能退化的初始时刻。

若V_c值曲线上升到某个时刻时后一时刻与该时刻之间的斜率达到最大值，则该时刻可判定为轴承失效时刻。此外，V_c值曲线在上升过程中的显著性转折点可认为是轴承性能退化过程中不同阶段的转折点。

与现有技术相比，本发明具有以下优点：

(1)本发明利用PCA对时域和频域组成的特征指标集进行加权融合，提取其第一主成分(PC1)作为滚动轴承性能退化评估的特征指标，主成分分析法(PCA)具有较好的空间映射和融合的特点，自身较为稳定，可以用来对时域和频域的特征指标进行加权融合，得到一个全面且有效的特征指标，能够用最小的特征信息表征主要的特征空间分布特性，从而能够更好的进行滚动轴承性能退化评估；

(2)本发明构造了基于SVDD算法的性能退化评估模型，相比于人工神经网络、隐马尔科夫等算法，本算法仅需要少量正常状态下的数据建模，克服了其他算法对故障模式下的数据的依赖，能够较好地解决故障诊断中的异常数据缺乏问题；

(3)本发明采用改进的HWPSO算法对SVDD模型的参数进行优化，既避免了标准鲸鱼算法易陷入局部最优的缺陷，又消除了人为选择模型参数的盲目性，从而最大限度地保证基于SVDD模型的轴承性能退化评估的准确性。

(4)本发明中通过在PSO算法中引入WOA算法进行空间搜索，充分了发挥了WOA算法的空间搜索能力，有效避免了PSO的“早熟”问题，同时利用PSO约束WOA的空间搜索范围，从而在更短的时间内找到全局最优值。

附图说明

图1为本发明实施例滚动轴承性能退化评估方法流程图；

图2为本发明实施例滚动轴承振动信号特征提取方法流程图；

图3为本发明实施例利用改进的鲸鱼粒子群混合算法得到最优参数的方法流程图；

图4为试验平台结构示意图；

图5为PCA约减融合后的特征指标；

图6为使用现有标准粒子群算法(PSO)和标准鲸鱼算法(WOA)寻优SVDD模型参数与本发明的结果对比图。

其中，1、轴承，2、电机，3、加速度传感器，4、热电偶。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。

实施例

如图1所示，本发明提供一种基于改进的HWPSO-SVDD的滚动轴承性能退化评估方法，该方法具体包括以下步骤：

步骤1：获取轴承正常状态下的振动信号x，并提取振动信号x的时域特征指标和频域特征指标。

其中，时域特征指标主要分为有量纲的时域特征指标和无量纲的时域特征指标，有量纲的时域特征指标与滚动轴承的实际工况相关联，包括滚动轴承的转速和载荷等；无量纲的时域特征指标不会受到滚动轴承不同工作条件的影响，它们具备直观反映滚动轴承性能退化趋势的能力。

时域特征指标可以在一定程度反映滚动轴承的性能退化趋势，但是这些特征指标的表现形式不同。

有量纲的时域特征指标包括均值和均方根值，通过均值可以看出振动信号的特征值是否稳定，如果特征值发生上下较大浮动，说明滚动轴承开始发生退化；均方根值又叫有效值，经常用于滚动轴承的故障检测和诊断，反映振动信号的幅值和能量，当幅值由平缓到不断发生上升的时候，说明滚动轴承开始退化。

无量纲的时域特征指标中的偏斜度和峭度指标的值会随着滚动轴承发生退化而变大；脉冲和裕度指标对早期的滚动轴承退化趋势较为敏感，幅值会发生明显上升的趋势，但是经过一段时间又会慢慢下降；峰值指标和偏斜度指标相对于其它的无量纲的特征指标，其稳定性更好，但是对于早期的滚动轴承的退化趋势不容易发现。所以这些时域特征指标在滚动轴退化过程中，都会有其相对应的表现形式，具体如表1和表2所示。

表1有量纲的时域特征指标

表2无量纲的时域特征指标

其中，x为振动信号，i＝1，2，3...N，N为振动信号的个数，x_i为第i个振动信号。

伴随着故障的发展，有量纲的时域特征指标会上升，而无量纲的时域特征指标则不会，它由概率密度函数来决定；在反映形式上，无论有量纲指标还是无量纲指标都是相对统一的，都能体现牵引电机滚动轴承的退化趋势，因此，用它们作为时域的特征指标。

通过振动信号的频谱，可以观察出振动能量随频率的变化状态，当滚动轴承发生退化趋势的时候，振动信号的频率成分的能量大小和主频位置都会改变，因此对于转向架牵引电机滚动轴承退化趋势的判断，频域分析相比时域分析，不仅仅只是局限于从滚动轴承的表面判断是否发生了退化，它还可以对滚动轴承运转时的特征故障频率进行分析，找出故障点，对滚动轴承的退化趋势进行预测。

对于原始的振动信号数据，提取其时域特征指标值，然后利用傅立叶变换将时域问题转换成频域问题。主要是将原来复杂的时间波形转换为相对简单直接的谐波来分析，容易得到信号的频率结构，可以计算出各谐波幅值、相位和能量与频率之间相互的数学关系式。表3为13个频域特征指标的数学表达式，其中p₁为中心频率；p₅为重心频率；p₄为频率方差；p₈为频域幅值平均值；p₁₀为均方根率：

表3频域特征指标

式中s(k)为振动信号x的频谱，k＝1，2，3...K，K是谱线数，f_k是第k条谱线数的频率值。p₁为中心频率，表示的是频域振动能量的大小；p₂-p₄、p₆和p₁₀-p₁₃表示的是频域集中或者分散程度；p₅和p₇-p₉表示的是主频带位置的变换。

步骤2：利用主成分分析(PCA)对时域特征指标与频域特征指标进行特征融合，得到表征滚动轴承运行趋势的融合特征指标。

主成分分析(PCA)是一种多元统计技术，它利用线性变换的方法将原始数据转化为少数几个主成分来表征原始数据的特征信息。PCA可以有效的处理原始数据空间维数过大的问题，通过加权融合将高维空间的数据降维至低维空间，使得原始数据的特征信息量变得简单和直观，因此PCA常被用来处理模式识别和图像处理等实际问题。

如图2所示，该步骤具体包括：

步骤21：构建振动信号的特征指标矩阵Y：

其中，y_ij为第j个特征指标中的第i个数据，n为每一种特征指标中的包含的数据个数，D为特征指标的种类数。

步骤22：计算特征指标矩阵Y中各特征指标数据的均值：

其中，

为第j个特征指标中数据均值。

步骤23：通过均值求解出方差

步骤24：对特征指标矩阵Y标准化，得到标准特征指标矩阵Y₀：

步骤25：计算协方差矩阵S：

步骤26：获取协方差矩阵S对角线上的各个特征值λ₁，λ₂，...，λ_n，并从大到小依次排列。各个特征值分别对应着各个主成分，表征振动信号x的特征信息，每个主成分对应特征值的其大小反应其所占据的总特征信息的比例，这个比例被称为贡献率，主成分的贡献率越大，说明其越能代表原始振动信号x的特征信息。

步骤27：计算各主成分的贡献率：

其中，其中为w_r第r个主成分对应的贡献率，λ_i为第i个特征值，L为所构造的主成分数。

步骤28：选择贡献率高于80％的主成分，作为融合特征指标。一般情况下第一主成分(PC1)便可以满足要求。

步骤3：将高斯核函数引入到支持向量数据描述算法(SVDD算法)中，利用HWPSO(Hybrid Whale-Particle Swarm Optimization Algorithm)优化算法对SVDD模型的惩罚参数C和核参数δ进行寻优。HWPSO优化算法是一种改进的基于种群混合的元启发算法，它结合了粒子群算法(PSO)和鲸鱼群算法(WOA)。

如图3所示，具体包括以下步骤：

步骤31：建立适应度函数F_sv：

其中，N_sv为SVDD模型里最小化支持向量个数，G为训练样本数。所以鲸鱼粒子群混合算法每尝试一次惩罚参数C和核参数δ，计算适应度，都要训练一次SVDD模型。

步骤32：初始化设置，具体为：

设置HWPSO最大迭代次数itermax、PSO粒子数pop、维数dim、解的上界和下界以及PSOgbest适应度和领头鲸适应度，其中PSOgbest适应度和领头鲸适应度均设为无穷大，并初始化单个粒子的局部最佳位置pbest，该值每次迭代都会更新。

步骤33：通过迭代，对SVDD模型的惩罚参数C和核参数δ寻优，迭代的过程具体包括：

步骤3301：计算惯性权重w、认知加速因子c₁、社会加速度因子c₂：

其中，w_min为惯性权重的最小值，其值为0.4，w_max为惯性权重的最大值，其值为0.9，itermax为HWPSO最大迭代次数，iter为当前HWPSO迭代次数，c_1i为认知加速度因子的初始值，c_1f为人质加速度因子的最终值，c_2i为社会加速度因子的初始值，c_2f为社会加速度因子的最终值，c_1i和c_2f的值为2.5，c_1f和c_2i的值为0.5。

步骤3302：由步骤31的公式计算每个粒子的当前适应度，并更新每个粒子的局部最佳位置pbest，若当前粒子适应度＜该粒子pbest适应度，则用当前适应度替换该粒子pbest适应度及其位置，若当前最优粒子pbest适应度＜PSOgbest适应度，则用该粒子pbest适应度替换PSOgbest适应度，粒子群每个粒子能够通过一定规则估计自身位置的适应值，每个粒子能够记住自己当前所找到的最好位置，称为局部最优pbest，对应PSOgbest适应度，此外还记住群体中所有粒子找到的一个最好位置，称为全局最优gbest，对应Gbest适应度。

步骤3303：更新每个粒子的速度和位置，其计算式为：

其中，

为粒子i在第k次迭代中第d维的当前速度，

为粒子i在第k次迭代中第d维的当前位置，

为粒子i在第k次迭代中第d维的更新后速度，即粒子i在第k+1次迭代中第d维的当前速度，

为粒子i在第k次迭代中第d维的更新后位置，即粒子i在第k+1次迭代中第d维的当前位置，rand₁为0到1之间随机数，pbest_id为粒子i第d维的局部最佳位置，rand₂为0到1之间随机数，gbest为全局最佳位置。

步骤3304：根据边界检查并修正新位置，该步骤中如果越界，则用最大值或最小值代替。

步骤3305：用粒子种群的最佳位置初始化鲸鱼种群，若gbest适应度＜领头鲸适应度，则为领头鲸指定gbest的位置和适应度。

步骤3306：计算WOA最大迭代次数im：

im＝[A×(iter)+B]

im为WOA最大迭代次数，A、B为常数，其值根据要优化的函数的模式来选择，本实施例中，

iter_max是HWPSO最大迭代次数，iter为HWPSO当前迭代次数，随着HWPSO的迭代，WOA的最大迭代次数im会逐渐减小，最后变为零。

步骤3307：获取每条鲸鱼的当前适应度，若当前适应度＜领头鲸适应度，则将该鲸鱼的当前适应度和位置指定为领头鲸的适应度和位置。

步骤3308：更新当前鲸鱼群个体的空间位置，具体过程为：

若概率因子p_m＜0.5，且系数向量A₂满足|A₂|≥1时，则更新鲸鱼群个体空间位置的表达式为：

X(t+1)＝X_rand(t)-A₂D₁

D₁＝|C₂X_rand(t)-X(t)|

若概率因子p_m＜0.5，且系数向量A₂满足|A₂|＜1时，则更新鲸鱼群个体空间位置的表达式为：

X(t+1)＝X^*(t)-A₂D₂

D₂＝|C₂X^*(t)-X(t)|

C₂＝2r

A₂＝2ar-a

其中，概率因子p_m为0到1的随机数，X(t+1)为更新的鲸鱼个体位置，X(t)为当前鲸群个体空间位置，X_rand(t)为当前鲸群中鲸鱼个体随机位置，X^*(t)为迄今鲸鱼个体最佳位置，t为当前WOA迭代次数，A₂和C₂为系数向量，D₁表示个体X(t)与随机位置X_rand(t)的距离，D₂表示包围步长，a为常数，并在迭代过程中从2减少到0，a_min＝0，a_max＝2，r为取值范围为[0，1]的随机向量。

若概率因子p_m≥0.5，则更新鲸鱼群个体空间位置的表达式为：

X(t+1)＝D₃e^blcos(2πl)+X^*(t)

D₃＝|X^*(t)-X(t)|

l＝(a₂-1)×rand+1

其中，b为常数，其值为1，l为介于[-1，1]之间的随机数，a₂为常数，D₃表示鲸群第i头鲸目前最佳位置到猎物之间的距离，rand为0到1间的随机数。

步骤3309：根据边界检查修正新位置，该步骤中如果越界，则用最大值或最小值代替，并判断当前WOA迭代次数t是否达到WOA最大迭代次数im，若未达到，则当前WOA迭代次数t的值加1，并返回执行步骤3307，否则进行下一步。

步骤3310：更新每个粒子的位置，若鲸鱼种群适应度＜PSO种群适应度，则将鲸鱼位置指定给粒子位置，若领头鲸适应度＜gbest适应度，则将领头鲸的位置指定给粒子种群最佳位置PSOgbest。

步骤3311：判断当前HWPSO迭代次数iter是否达到HWPSO最大迭代次数iter_max，若是，则执行下一步，否则当前HWPSO迭代次数iter的值加1，并返回执行步骤3301。

步骤3312：输出粒子种群最佳位置PSOgbest，其对应的坐标即是SVDD模型最优的惩罚参数C和核参数δ。

步骤4：根据寻优完成的惩罚参数C和核参数δ，构建优化的SVDD模型，即HWPSO-SVDD模型，并将步骤2得到的融合特征指标作为训练样本，输入经过HWPSO参数优化的SVDD模型进行模型训练，得到训练完成的性能退化评估模型，并求出超球体的半径R。

SVDD评估模型主要是通过输入特征指标建立SVDD模型，SVDD的基本思想就是生成一个最小超球体，使其尽可能包含所有的正常特征样本。优化目标为寻找在特征空间包含所有样本点的半径最小的超球面，其目标函数表示为：

s.t.(z_i-c)^T(z_i-c)≤R²+ε_i，ε_i≥0，i＝1，2，…，M

其中，z_i为训练样本，M为训练样本个数，C为惩罚因子，ε_i为松弛变量，c为超球面中心，R为超球面半径，K(·)为核函数，α_i为拉格朗日乘子，z_sv为支持向量(位于超球面上的样本称为支持向量)。

步骤5：对待测的轴承振动数据根据步骤2的方法得到待测融合特征指标，并输入性能退化评估模型，计算各融合特征指标距离超球体中心的距离，得到性能退化指标，即V_c值，V_c值的计算式为：

其中，V_c为融合特征指标距离超球体中心的距离，z_q为待测样本。

步骤6：设置自适应报警阈值线，根据V_c值对轴承的早期退化状态做出及时准确的预警。

V_c是一个连续变化的参数，表示设备偏离正常状态的程度。设定报警阈值，能够对轴承的健康状态进行实时的监控。根据统计学中的3σ法则，对于一个均值为

方差为σ²的高斯随机变量而言，取值落在区间

内的概率为99.73％。一旦某个取值超过该范围，有理由认为该取值来自其他变量。同样可以假设处于相近性能退化状态下的V_c值也符合一定的正态分布，一旦连续多个V_c值超过前面的V_c值界定的3σ取值范围(单个V_c值的超限可能是受外界干扰影响)，则认为轴承的性能退化状态已发生较大改变。该阈值不需要事先设定，根据实际所测得的V_c值动态改变，是自适应的。由于V_c是性能退化程度的增函数，所以只考虑阈值上限，即

自适应报警阈值的计算公式为：

其中，t_s表示可以确定的轴承处于正常状态的时刻，T(t)表示t时刻的V_c值，mean、std分别表示求均值和标准偏差。阈值T_h的求解分为3个阶段：第1阶段数据来源于早期无故障状态，是一个固定值；第2阶段将t时刻的T(t)与t-1时刻的T_h(t-1)进行比较，若T(t)在T_h(t-1)范围内，则将T(t)纳入原始数据，计算T_h(t)；第3阶段为其后连续N_u个V_c值均超限，则定义t＝t_e为轴承性能退化的初始时刻。

若V_c值曲线上升到某个时刻时后一时刻与该时刻之间的斜率达到最大值，则该时刻可判定为轴承失效时刻。此外，V_c值曲线在上升过程中的显著性转折点可认为是轴承性能退化过程中不同阶段的转折点。

实施例1

本实施例中，采用来源于美国辛辛那提大学智能维护中心的试验数据进行试验，试验平台如图4所示。

电机2带动主轴以2000r/min转速运动，轴承试验台与电机2连接的转轴上安装了四个Rexnord ZA-2115滚动轴承，滚子直径为8.407mm，节圆直径为71.501mm，接触角为15.17°，每个轴承1的径向方向上安装一个加速度传感器3，并且安装有热电偶4，在轴和轴承上加载约26.67kN的径向载荷，四个轴承1均采用油润滑。轴承振动信号由NIDAQ 6062E数据采集卡每隔10min采集一次。采样频率为20kHz，每个数据文件的采样长度为20480个点。在运行了约163h后，轴承2因出现严重外圈故障而失效，期间共采集了984组数据，此数据集即为轴承2的全寿命周期试验数据。本节采用轴承2的全生命周期数据进行性能退化评估试验。

由于时域、频域等包含的指标量较多，采用单个指标量无法全面地反映滚动轴承的退化趋势，且单个时域或频域指标还存在评估能力不足的问题。因此，为提取一个能有效地表征滚动轴承退化趋势的特征量，利用PCA融合时域指标和频域指标。提取其第一主成分(PC1)作为滚动轴承的退化趋势性能预测指标代入HWPSO-SVDD模型中，训练SVDD模型并完成滚动轴承的性能退化评估。利用PCA融合时域指标和频域指标，提取其第一主成分如图5所示。由图5可知，选PCA第一主成分(PC1)作为轴承性能退化特征指标时，在500点到700点之间开始轻微上升，700之后有着明显的变化趋势，所以选用500之前的数据进行训练。

进行滚动轴承性能退化评估时，首先利用前300组正常状态下的数据作为正常样本，训练SVDD模型，利用改进的HWPSO算法优化惩罚参数和核参数，取种群规模N_size＝10，最大迭代次数g_max＝10，最终寻找到的最优惩罚参数C和核参数δ。

得到SVDD超球体后，将984组数据全部作为待测样本，依照上述方法输入到SVDD模型中，得到轴承全生命周期的V_c值。为了减少外界干扰对V_c值的影响，对其进行五点平滑处理。如果连续有7个时刻的V_c值超限，则认为轴承已产生了早期故障，即N_u＝7。为验证所提方法的优越性，采用基于PSO-SVDD和WOA-SVDD的性能退化评估方法对轴承2的全寿命数据进行评估，评估结果如图6所示。从图6中可以看出，最早故障出现在第634个和第633个样本处，与HWPSO算法优化的模型相比晚了370分钟左右(每隔10min采集一次数据)，且在最早故障出现后曲线的上升趋势并没有比HWPSO-SVDD模型明显。通过对比可见，本发明所提方法在早期故障检测方面及性能退化评估方面都具有优越性，下表为3种方法预测结果与真实值之间的误差比较结果：

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的工作人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (10)

1.一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，包括以下步骤：

S1：获取正常滚动轴承振动信号的特征指标，通过PCA算法得到融合特征指标，建立训练样本；

S2：利用HWPSO优化算法对SVDD模型的惩罚参数C和核参数δ进行寻优；

S3：利用寻优得到的惩罚参数C和核参数δ，构建HWPSO-SVDD模型；

S4：将训练样本输入HWPSO-SVDD模型进行模型训练，得到训练完成的性能退化评估模型；

S5：获取待测滚动轴承振动信号的融合特征指标，利用训练完成的性能退化评估模型，计算得到待测滚动轴承的性能退化指标V_c值；

S6：设置自适应报警阈值，当V_c值超过自适应报警阈值时，对轴承的早期退化状态做出预警。

2.根据权利要求1所述的一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，所述的滚动轴承振动信号的特征指标包括时域特征指标和频域特征指标，所述的时域特征指标包括有量纲时域特征指标和无量纲时域特征指标，所述的频域特征指标包括中心频率、重心频率、频率方差、频域幅值平均值和均方根率。

3.根据权利要求2所述的一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，所述的步骤S1具体包括：

S11：构建特征指标矩阵Y：

其中，y_ij为第j个特征指标中的第i个数据，n为每一种特征指标中的包含的数据个数，D为特征指标的种类数；

S12：计算特征指标矩阵Y中各特征指标的均值：

其中，

为第j个特征指标中数据均值；

S13：求解特征指标的方差

S14：对特征指标矩阵Y标准化，得到标准特征指标矩阵Y₀：

S15：计算特征指标协方差矩阵S：

S16：获取特征指标协方差矩阵S对角线上的各特征值λ₁,λ₂,...,λ_n，并从大到小依次排列；

S17：计算各主成分的贡献率：

其中，其中为w_r第r个主成分对应的贡献率，λ_i为第i个特征值，L为所构造的主成分数；

S18：选择贡献率高于80％的主成分，作为融合特征指标。

4.根据权利要求1所述的一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，所述的步骤S2具体包括：

S21：建立适应度函数F_sv；

S22：初始化设置改进粒子群算法参数和鲸鱼群算法参数，并分别设置PSOgbest适应度和领头鲸适应度为无穷大；

S23：通过改进粒子群算法和鲸鱼群算法的迭代，对SVDD模型的惩罚参数C和核参数δ寻优，得到最优解。

5.根据权利要求4所述的一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，所述的适应度函数F_sv的表达式为：

其中，N_sv为SVDD模型里最小化支持向量个数，G为训练样本数。

6.根据权利要求4所述的一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，所述的步骤S23具体包括：

S231：计算惯性权重w、认知加速因子c₁、社会加速度因子c₂：

其中，w_min为惯性权重的最小值，其值为0.4，w_max为惯性权重的最大值，其值为0.9，iter_max为HWPSO最大迭代次数，iter为当前HWPSO迭代次数，c_1i为认知加速度因子的初始值，c_1f为认知加速度因子的最终值，c_2i为社会加速度因子的初始值，c_2f为社会加速度因子的最终值，c_1i和c_2f的值为2.5，c_1f和c_2i的值为0.5；

S232：获取每个粒子的当前适应度，判断粒子当前适应度是否小于该粒子pbest适应度，若是，则用该粒子当前适应度及其位置替换该粒子pbest适应度及其位置；

S233：判断当前最优粒子pbest适应度是否小于PSOgbest适应度，若是，则用该粒子pbest适应度替换PSOgbest适应度；

S234：更新每个粒子的速度和位置；

S235：根据边界检查修正更新全局最佳位置，若修正的最新全局最佳位置越界，则用最大值或最小值代替；

S236：利用粒子种群的最佳位置初始化鲸鱼种群算法，并通过鲸鱼种群算法的迭代，对粒子种群最佳位置PSO gbest进行优化；

S237：判断当前HWPSO迭代次数iter是否达到HWPSO最大迭代次数iter_max，若是，则执行步骤S238，否则当前HWPSO迭代次数iter的值加1，并返回执行步骤S231；

S238：输出粒子种群最佳位置PSO gbest，其对应的坐标即是SVDD模型最优的惩罚参数C和核参数δ。

7.根据权利要求6所述的一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，所述的步骤S236具体包括：

S2361：判断gbest适应度是否小于领头鲸适应度，若是，则将gbest适应度及其位置设定为领头鲸适应度和位置；

S2362：计算WOA最大迭代次数im：

im＝[A×(iter)+B]

其中，A、B为常数；

S2363：获取每条鲸鱼的当前适应度，判断鲸鱼当前适应度是否小于领头鲸适应度，若是，则将该鲸鱼当前适应度及其位置设定为领头鲸的适应度和位置；

S2364：更新当前鲸鱼群个体的空间位置；

S2365：根据边界检查修正更新领头鲸位置，若修正的领头鲸位置越界，则用最大值或最小值代替；

S2366：判断当前WOA迭代次数t是否达到WOA最大迭代次数im，若是，则执行步骤S237，否则当前WOA迭代次数t的值加1，并返回执行步骤S2363。

8.根据权利要求7所述的一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，所述的更新每个粒子的速度和位置的计算式为：

其中，

为粒子i在第k次迭代中第d维的当前速度，

为粒子i在第k次迭代中第d维的当前位置，

为粒子i在第k次迭代中第d维的更新后速度，即粒子i在第k+1次迭代中第d维的当前速度，

为粒子i在第k次迭代中第d维的更新后位置，即粒子i在第k+1次迭代中第d维的当前位置，rand₁为0到1之间随机数，pbest_id为粒子i第d维的局部最佳位置，rand₂为0到1之间随机数，gbest为全局最佳位置；

所述的更新当前鲸鱼群个体的空间位置，具体包括：

若概率因子p_m＜0.5，且系数向量A₂满足|A₂|≥1时，则更新鲸鱼群个体空间位置的表达式为：

X(t+1)＝X_rand(t)-A₂D₁

D₁＝|C₂X_rand(t)-X(t)|

若概率因子p_m＜0.5，且系数向量A₂满足|A₂|＜1时，则更新鲸鱼群个体空间位置的表达式为：

X(t+1)＝X^*(t)-A₂D₂

D₂＝|C₂X^*(t)-X(t)|

C₂＝2r

A₂＝2ar-a

其中，概率因子p_m为0到1的随机数，X(t+1)为更新的鲸鱼个体位置，X(t)为当前鲸群个体空间位置，X_rand(t)为当前鲸群中鲸鱼个体随机位置，X*(t)为迄今鲸鱼个体最佳位置，t为当前WOA迭代次数，A₂和C₂为系数向量，D₁表示个体X(t)与随机位置X_rand(t)的距离，D₂表示包围步长，a为常数，并在迭代过程中从2减少到0，a_min＝0，a_max＝2，r为取值范围为[0,1]的随机向量；

若概率因子p_m≥0.5，则更新鲸鱼群个体空间位置的表达式为：

X(t+1)＝D₃e^blcos(2πl)+X^*(t)

D₃＝|X^*(t)-X(t)|

l＝(a₂-1)×rand+1

其中，b为常数，其值为1，l为介于[-1，1]之间的随机数，a₂为常数，D₃表示鲸群第i头鲸目前最佳位置到猎物之间的距离，rand为0到1间的随机数。

9.根据权利要求1所述的一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，所述的HWPSO-SVDD模型的目标函数为：

s.t.(z_i-c)^T(z_i-c)≤R²+ε_i,ε_i≥0,i＝1,2,...,M

其中，z_i为目标样本，M为训练样本个数，C为惩罚因子，ε_i为松弛变量，c为超球面中心，R为超球面半径，K(·)为核函数，α_i为拉格朗日乘子，z_sv为支持向量。

10.根据权利要求9所述的一种基于HWPSO-SVDD模型的滚动轴承性能退化评估方法，其特征在于，所述的待测滚动轴承的性能退化指标V_c值通过计算待测数据特征样本q距离超球中心的距离得到，其计算式为：

其中，V_c为融合特征指标距离超球体中心的距离，z_q为待测样本；

所述的自适应报警阈值T_h(t)根据实际所测得的V_c值自适应改变，其计算公式为：

其中，t_s为确定的轴承处于正常状态的时刻，T(t)为t时刻的V_c值，mean、std分别表示求均值和标准偏差，t_e为轴承性能退化的初始时刻。

Images