CN106644481A - 基于数学形态学和ifoa‑svr的滚动轴承可靠度预测方法 - Google Patents

Abstract

基于数学形态学和IFOA‑SVR的滚动轴承可靠度预测方法，涉及滚动轴承可靠度预测技术领域。为保证预测精度同时增加预测步长而提出的预测方法。该方法首先提取振动信号的包络信号，计算该包络信号的数学形态学分形维数，将其作为滚动轴承性能退化状态特征；其次，利用IFOA对SVR中的参数C，g以及ε同时进行寻优，建立预测模型。同时，利用极大似然估计结合IFOA建立威布尔比例故障率模型，进而得到可靠度模型；最后，将退化状态特征作为IFOA‑SVR预测模型的输入，采用长期迭代预测法获取特征预测结果，并将该结果嵌入到可靠度模型中，从而预测出轴承运行状态的可靠度。实验表明，利用所提方法在保证预测精度的前提下增加预测步长。

Description

基于数学形态学和IFOA-SVR的滚动轴承可靠度预测方法

技术领域

本发明涉及一种基于数学形态学和IFOA-SVR的滚动轴承可靠度预测方法，涉及滚动轴承可靠度预测技术领域。

背景技术

滚动轴承是旋转机械中的关键部件，一旦出现故障将会造成大量的经济损失甚至危害人的生命安全^[1-2]。因此，准确预测滚动轴承在下一阶段的工作状态是合理制定机械设备维修计划的前提基础^[3-4]。

目前，滚动轴承振动信号的特征提取方法研究受到学者的广泛关注。文献[5]提出了基于形态分量分析和包络谱的滚动轴承故障诊断方法，该方法可有效提取滚动轴承振动信号中的故障特征。文献[6]提出局部均值分解与形态学分形维数相结合的滚动轴承故障诊断方法，可有效地对滚动轴承进行故障诊断。本发明将Hilbert包络解调与形态学分形维数结合，对滚动轴承全寿命数据提取出数学形态学分形维数，并将其作为性能退化状态特征。

在建立支持向量回归(Support vector regression,SVR)模型方面，文献[7]利用粒子群算法(Particle swarm optimization,PSO)优化最小二乘SVR对滚动轴承的性能退化趋势进行了20点预测，平均绝对百分误差(Mean absolute percentage error,MAPE)为15.82％。文献[8]利用交叉验证法对C，g寻优建立SVR模型，对滚动轴承故障趋势进行了7点预测，MAPE为0.2％。文献[9]通过经验选取C，g参数建立SVR模型，并对电厂汽轮发电机组转子的振动幅值进行了30点预测，平均绝对误差(Mean absolute error,MAE)为0.6424，MAPE为1.48％。文献[10]利用PSO对SVR模型中的C，g以及ε寻优，预测航空飞行器的剩余寿命。文献[7]到文献[9]预测步长较短，文献[8]和[9]只对SVR模型中的C和g进行选取，并没有提出对不敏感误差ε的取值问题。同时，凭经验选取C和g有很大的盲目性，通过交叉对比测试确定C和g的值则需要进行大量的实验，而PSO相对复杂，并且容易陷入局部最优，导致对模型参数的赋值不够准确。文献[11]提出一种果蝇优化算法(Fruit fly optimizationalgorithm,FOA)。文献[12]提出递减步长果蝇优化算法(Diminishing step fruit flyoptimization algorithm,DS-FOA)实现了全局搜索能力和局部寻优能力的平衡。但DS-FOA在二维空间进行搜索，不能真实反映果蝇的觅食行为。SVR模型内部参数的选取会对SVR的预测精度以及预测步长产生较大影响，本发明采用在三维空间搜索的改进果蝇优化算法(Improved fruit fly optimization algorithm,IFOA)对SVR模型中的3个参数C、g以及ε同时进行参数寻优。并且在进行预测时，将MAE、均方根误差(Root mean square error,RMSE)、归一化均方误差(Normalized mean square error,NMSE)以及MAPE的和作为适应度函数，相比于任一误差作为适应度函数具有更强的寻优能力。

在可靠度模型方面，可靠度理论已被应用于机械工程以及航空航天等各个领域^[13]。比例故障率模型的优势在于将轴承的性能退化特征与可靠度理论相结合，这样，轴承的可靠度评估就能够在设备状态监测数据的基础上得到更新。文献[14]提出将峭度、均方根作为特征指标，利用威布尔比例故障率模型(Weibull proportional hazard model,WPHM)作为可靠度评估模型，并利用fminsearch优化函数求解极大似然方程组确定WPHM的待定参数，有效地对铁路机车轮对轴承进行可靠度评估。文献[15]在估计WPHM的待定参数时，利用牛顿迭代法对极大似然方程组进行求解。但上述两种求解方法都需要根据经验设定搜索初值，并且运算时间较长。本发明将极大似然方程组中每个方程的绝对值的和作为IFOA的适应度函数，向适应度函数最小的方向优化，可以快速求出方程组的解，确定WPHM的待定参数。

发明内容

本发明为了保证滚动轴承运行状态可靠度的预测精度同时增加预测步长，提出一种数学形态学分形维数结合改进果蝇优化算法-支持向量回归(Improved fruit flyoptimization algorithm-support vector regression,IFOA-SVR)的滚动轴承可靠度预测方法，即提出一种基于数学形态学和IFOA-SVR的滚动轴承可靠度预测方法。

本发明为解决上述技术问题采取的技术方案是：

一种基于数学形态学和IFOA-SVR的滚动轴承可靠度预测方法，所述方法的实现过程为：

步骤一、求出包络信号的数学形态学分形维数D_M；

步骤二、从获取的数学形态学分形维数D_M选取训练样本对，基于所述训练样本对利用IFOA对SVR模型中的参数C、g以及ε同时进行寻优，建立IFOA-SVR预测模型：

构建IFOA的过程：

(1)初始化算法参数：设置果蝇种群规模Sizepop，最大觅食代数Maxgen，并随机初始化果蝇群体位置坐标(X₀,Y₀,Z₀)；

(2)果蝇个体利用嗅觉随机搜索的方向和距离可以通过式(12)获得

式中，i＝1,2,...,Sizepop，L₀为初始步长值，gen为当前觅食代数；

(3)由于无法确定食物源的具体位置，所以需要通过式(13)估计第i个果蝇个体的当前位置与坐标原点间的距离Dist_i，之后计算出味道浓度判定值S_i

S_i＝1/Dist_i (14)

(4)将S_i代入味道浓度判定函数，计算出果蝇个体当前位置的味道浓度

Smell_i＝function(S_i) (15)

(5)当前果蝇群体中具有最高味道浓度的个体，可由式(16)获得

[bestSmell,bestIndex]＝max(Smell_i) (16)

式中，bestSmell表示果蝇群体中具有最高味道浓度的个体的味道浓度值，bestIndex表示果蝇群体中具有最高味道浓度的个体的位置；

(6)保留果蝇群体中最佳味道浓度值和与其对应的个体坐标，同时果蝇群体利用自身的视觉对食物源进行定位，然后飞往食物源所在的位置；

(7)进入迭代寻优过程，重复步骤(2)-(5)，并判断当前味最高道浓度是否好于前一迭代味道浓度，且gen<Maxgen；若成立，则执行步骤(6)；

建立IFOA-SVR模型，其过程为：

(1)初始化IFOA参数，包括Sizepop、Maxgen以及果蝇个体初始位置；对SVR中的3个参数：惩罚系数C、核函数参数g以及不敏感误差ε进行寻优，所以初始坐标为(X₀ ¹,Y₀ ¹,Z₀ ¹)，(X₀ ²,Y₀ ²,Z₀ ²)和(X₀ ³,Y₀ ³,Z₀ ³)；

(2)附与每个果蝇个体随机方向和飞行距离，并用递减搜索步长代替固定步长，得到(X_i ¹,Y_i ¹,Z_i ¹)，(X_i ²,Y_i ²,Z_i ²)，(X_i ³,Y_i ³,Z_i ³)，计算当前果蝇个体与原点间距离的倒数，得到味道浓度判定值S_i ¹，S_i ²和S_i ³；

(3)确定SVR中参数C、g和ε的范围，即C∈[2^-14,2¹⁴]，g∈[2^-14,2¹⁴]，ε∈[2^-14,2¹⁴]；

(4)将所述训练样本对输入到SVR模型中，对模型进行训练，将平均绝对误差MAE、均方根误差RMSE、归一化均方误差NMSE以及平均绝对百分误差MAPE的和作为适应度函数，即

(5)找到适应度函数Fitness的最小值对应的果蝇个体，开始迭代寻优，并判断最小Fitness是否低于前一代最小Fitness；如低于，则保留最小Fitness值及其对应的坐标，并将其赋给初始坐标；如高于，则返回步骤(2)；

(6)找到C、g和ε的最佳值，建立IFOA-SVR预测模型；

步骤三、从所述性能退化状态特征中截取一定长度Z作为WPHM的响应协变量，结合极大似然估计获得似然函数方程组，把方程组中每个方程绝对值的和作为IFOA的适应度函数，求出方程组的解，确定WPHM的待定参数，进而得到可靠度模型R(t,Z)，过程如下：

基于WPHM建立滚动轴承性能退化状态特征与可靠度之间的数学关系，WPHM的表达式为

式中β为形状参数，η为尺度参数，μ为协变量回归参数，t表示时间，Z为截取一定长度的性能退化状态特征，即从数学形态学分形维数D_M中截取一定长度数据；

h(t,Z)和可靠度函数R(t,Z)之间的关系为

从而，R(t,Z)可表示为

利用极大似然估计法得到关于待定参数β、η、μ的方程组并利用IFOA对该方程组进行求解，确定WPHM的表达式h(t,Z)的待定参数β、η、μ的最终取值；将参数β、η、μ代入式(21)即可建立可靠度模型；

步骤四、将性能退化状态特征的最后一部分数据作为IFOA-SVR预测模型的输入，采用长期迭代预测法获取性能退化状态特征预测结果；

步骤五、将IFOA-SVR预测模型得出的预测结果作为响应协变量嵌入到可靠度模型R(t,Z)中，即可计算出所述性能退化状态特征预测结果所对应的可靠度，实现滚动轴承可靠度的预测。

在步骤一中，通过Hilbert变换计算出状态监测数据的包络信号；

步骤一一、获取滚动轴承的状态监测数据：滚动轴承的振动信号；

步骤一二、通过Hilbert变换计算出状态监测数据的包络信号；

步骤一三、利用式(4)求出每次采集振动信号的包络信号的数学形态学分形维数D_M，将D_M作为滚动轴承性能退化状态特征；

式中：D_M为信号的Minkowski-Bouligand维数即数学形态学分形维数，c为常数，对和log(1/λ)进行最小二乘线性拟合即可得到对信号数学形态学分形维数D_M；λ＝1,2,…,λ_max，λ_max为最大分析尺度；

A_Se(λ)为定义在分析尺度λ下f(n)关于Se(n)膨胀和腐蚀运算的覆盖面积，f(n)表示滚动轴承振动信号的包络信号；Se(n)为结构元素，为一维离散向量；Se(n)的长度及λ_max的值通过控制变量法求得。

本发明的有益效果是：

本发明提出的一种基于数学形态学分形维数结合改进果蝇算法-支持向量回归(Improved fruit fly optimization algorithm-support vector regression,IFOA-SVR)的滚动轴承可靠度预测方法，充分发挥IFOA-SVR预测模型的优势，对提取的滚动轴承状态特征数据进行预测。将预测结果嵌入到已建立的可靠度模型中进行滚动轴承可靠度预测，在保证滚动轴承运行状态可靠度预测精度的同时增加了预测步长。

该方法首先提取振动信号的包络信号，计算该包络信号的数学形态学分形维数，将其作为滚动轴承性能退化状态特征；其次，利用IFOA对SVR中的参数C，g以及ε同时进行寻优，建立IFOA-SVR预测模型。同时，利用极大似然估计结合IFOA建立威布尔比例故障率模型(Weibull proportional hazard model,WPHM)，进而得到可靠度模型；最后，将退化状态特征作为IFOA-SVR预测模型的输入，采用长期迭代预测法获取特征预测结果，并将该结果嵌入到可靠度模型中，从而预测出轴承运行状态的可靠度。实验结果表明，利用所提方法对滚动轴承可靠度进行预测，能在保证预测精度的前提下增加预测步长。

具体而言，本发明将Hilbert包络解调与形态学分形维数结合，对滚动轴承全寿命数据提取出数学形态学分形维数，并将其作为性能退化状态特征。本发明采用在三维空间搜索的改进果蝇优化算法(Improved fruit fly optimization algorithm,IFOA)对SVR模型中的3个参数C、g以及ε同时进行参数寻优。并且在进行预测时，将MAE、均方根误差(Rootmean square error,RMSE)、归一化均方误差(Normalized mean square error,NMSE)以及MAPE的和作为适应度函数，相比于任一误差作为适应度函数具有更强的寻优能力。本发明将极大似然方程组中每个方程的绝对值的和作为IFOA的适应度函数，向适应度函数最小的方向优化，可以快速求出方程组的解，确定WPHM的待定参数。

附图说明

图1是本发明中的果蝇三维空间觅食示意图，图2是本发明的可靠度预测方法流程图；

图3为本发明方法在应用中的第860次采集数据的时域波形图，图4是本发明方法在应用中对应的第860次采集数据的包络信号波形图，图5为本发明在应用中第1点至第860点的可靠度曲线图，图6为本发明在应用中GA优化SVR模型得出的预测结果图，图7为本发明在应用中PSO优化SVR模型得出的预测结果图，图8为本发明在应用中DS-FOA优化SVR模型得出的预测结果图，图9为本发明在应用中IFOA优化SVR模型得出的预测结果图，图10为本发明在应用中50点预测数据的可靠度趋势曲线图；

图11为由包络求分形维数与直接求分形维数的对比图，本发明采用由由包络求分形维数。

具体实施方式

具体实施方式一：如图1至11所示，本实施方式针对所述的基于数学形态学和IFOA-SVR的滚动轴承可靠度预测方法的实现过程具体描述如下：

1基于数学形态学的分形维数

数学形态学包括两种基本算子，分别是膨胀运算和腐蚀运算。设原始信号f(n)和结构元素Se(n)分别为定义在集合F＝{0,1,…,N-1}和集合G＝{0,1,…,M-1}上的两个一维离散函数，且N≥M。在每一个分析尺度λ下，令Se(n)对f(n)进行一次膨胀和腐蚀运算，即：

式中：表示膨胀运算，Θ表示腐蚀运算，λ＝1,2,…,λ_max，λ_max为最大分析尺度。

定义在尺度λ下f(n)关于Se(n)膨胀和腐蚀运算的覆盖面积A_Se(λ)为

Maragos证明，A_Se(λ)满足式(4)

式中D_M为信号的Minkowski-Bouligand维数即数学形态学分形维数，c为常数，对和log(1/λ)进行最小二乘线性拟合即可得到对信号数学形态学分形维数的估计。

2支持向量回归机

SVR算法一般用于解决数据回归预测问题，通过建立一个最优超平面，使状态空间内的各个数据点距离该超平面最近，并将该超平面作为回归模型进行预测。

训练样本对(x₁,y₁)，(x₂,y₂)，…，(x_r,y_r)，x_i∈Rⁿ，y_i∈R，x_i为输入样本，y_i为输出样本，i＝1,2,...,r。SVR的回归函数表达式为

f(x)＝<ω·x_i>+b (5)

式中<ω·x_i>为ω和x_i的内积。系数ω和b通过求解式(6)的最小值优化问题获得

由于拟合曲线必然会存在误差，但要把误差控制在一定的可容许范围内，ε为不敏感误差，式(6)中L_ε为不敏感损失函数，其表达式为

对于回归错误的数据点，需要引入松弛变量ξ和ξ^*。将不敏感损失函数Lε代入式(6)可得

满足约束条件

求解上式时，一般采用对偶理论将其转化为二次规划问题。对于非线性数据，引入非线性映射函数Φ，建立拉格朗日方程，通过化简，可得式(8)的对偶式

式中α，α^*为拉格朗日乘子。

令K(x_i,x_j)＝Φ(x_i)^TΦ(x_j)为特征空间的内积，称K(x_i,x_j)为核函数，选择径向基函数作为核函数。根据KKT定理，可求得变量α_i，α_i ^*，b，最终可得支持向量机的回归函数为

3基于IFOA的SVR寻优

3.1 IFOA算法

FOA是一种基于果蝇觅食行为推演出寻求全局优化的算法。利用在三维空间搜索的IFOA对SVR模型进行优化，可以有效地获取模型中的最优参数。

如图1所示为果蝇三维空间觅食示意图。

IFOA步骤：

(1)初始化算法参数。设置果蝇种群规模Sizepop，最大觅食代数Maxgen，并随机初始化果蝇群体位置坐标(X₀,Y₀,Z₀)。

(2)果蝇个体利用嗅觉随机搜索的方向和距离可以通过式(12)获得

式中，i＝1,2,...,Sizepop，L₀为初始步长值，gen为当前觅食代数。

(3)由于无法确定食物源的具体位置，所以需要通过式(13)估计第i个果蝇个体的当前位置与坐标原点间的距离Dist_i，之后计算出味道浓度判定值S_i

S_i＝1/Dist_i (14)

(4)将S_i代入味道浓度判定函数，计算出果蝇个体当前位置的味道浓度

Smell_i＝function(S_i) (15)

(5)当前果蝇群体中具有最高味道浓度的个体，可以由式(16)获得

[bestSmell,bestIndex]＝max(Smell_i) (16)

(6)保留果蝇群体中最佳味道浓度值和与其对应的个体坐标，同时果蝇群体利用自身的视觉对食物源进行定位，然后飞往食物源所在的位置。

(7)进入迭代寻优过程，重复步骤(2)-(5)，并判断当前味最高道浓度是否好于前一迭代味道浓度，且gen<Maxgen；若成立，则执行步骤(6)。

3.2 IFOA-SVR模型

(1)初始化IFOA参数。包括Sizepop、Maxgen以及果蝇个体初始位置。由于本发明需要对SVR中的3个参数C、g以及ε进行寻优，所以初始坐标为(X₀ ¹,Y₀ ¹,Z₀ ¹)，(X₀ ²,Y₀ ²,Z₀ ²)和(X₀ ³,Y₀ ³,Z₀ ³)。

(2)附与每个果蝇个体随机方向和飞行距离，并用递减搜索步长代替固定步长，得到(X_i ¹,Y_i ¹,Z_i ¹)，(X_i ²,Y_i ²,Z_i ²)，(X_i ³,Y_i ³,Z_i ³)。计算当前果蝇个体与原点间距离的倒数，得到味道浓度判定值S_i ¹，S_i ²和S_i ³。

(3)确定SVR中参数C、g和ε的范围，即C∈[2^-14,2¹⁴]，g∈[2^-14,2¹⁴]，ε∈[2^-14,2¹⁴]。

(4)将训练样本对输入到SVR模型中，对模型进行训练，将MAE、RMSE、NMSE以及MAPE的和作为适应度函数，即

(5)找到Fitness的最小值对应的果蝇个体，开始迭代寻优，并判断最小Fitness是否低于前一代最小Fitness。如低于，则保留最小Fitness值及其对应的坐标，并将其赋给初始坐标。如高于，则返回步骤(2)。

(6)找到C、g和ε的最佳值，建立IFOA-SVR预测模型。

4滚动轴承运行状态可靠度

WPHM建立了设备运行状态特征指标与可靠度之间的数学关系，表达式为

式中β为形状参数，η为尺度参数，μ为协变量回归参数，Z为振动特征指标数据。

h(t,Z)和可靠度函数R(t,Z)之间的关系为

从而，R(t,Z)可表示为

利用极大似然估计法对待定参数β、η、μ进行估计，确定待定参数最终取值。

5滚动轴承可靠度预测方法及流程

在WPHM参数确定的情况下，建立可靠度模型，将特征预测结果嵌入已建立的可靠度模型中，即可预测未来一段时间内轴承运行的可靠度，具体的预测流程见图2。

(1)通过Hilbert变换计算出状态监测数据的包络信号，利用式(4)，求出包络信号的数学形态学分形维数D_M。

(2)选取适当的训练样本对，并利用IFOA对SVR模型中的参数C、g以及ε同时进行寻优，建立IFOA-SVR预测模型。

(3)选取一定长度的振动特征指标数据作为WPHM的响应协变量，结合极大似然估计获得似然函数方程组，把方程组中每个方程绝对值的和作为IFOA的适应度函数，求出方程组的解，确定WPHM的待定参数，进而得到可靠度模型。

(4)将振动特征指标数据作为IFOA-SVR预测模型的输入，采用长期迭代预测法^[16]获取振动特征预测结果。

(5)将IFOA-SVR预测模型得出的预测结果作为响应协变量嵌入到可靠度模型中，即可计算出这些特征预测结果所对应的可靠度，实现滚动轴承可靠度的预测。

针对本明提出的方法进行具体应用与分析

滚动轴承全寿命数据来自于Cincinnati大学^[17]，每隔10min进行一次数据采集，共采集984次，每次采样长度为20480点，采样频率为20KHz。图3为第860次采集数据的时域波形，图4为其对应的包络信号波形。

在每次采集的数据中，取中间的4096点进行Hilbert变换，求出包络信号，进而计算出数学形态学分形维数，并作为一个特征数据点。通过控制变量法实验对比证明，当λ_max＝11，Se(n)＝(0,0,0,0,0,0,0,0,0,0,0,0)时，利用式(4)求出包络信号的数学形态学分形维数能够作为刻画滚动轴承性能退化过程的特征。

将滚动轴承性能退化状态特征作为响应协变量，取600点到860点之间的特征指标值，采用极大似然估计，利用IFOA求得WPHM模型的待定参数，如表1所示。

表1 WPHM待定参数的估计结果

将待定参数代入式(19)，即可确定当数学形态学分形维数作为响应协变量时的WPHM，进而建立可靠度模型，按照式(22)可以计算出滚动轴承在任意时刻下的可靠度

式中m＝1,2,…,984，984表示数据采集总次数；

每个数据点等于10分钟，数据点与时间存在对应关系。

将第1点至第860点之间对应的特征指标值嵌入可靠度模型即可计算出在这段时间内轴承运行的可靠度，如图5所示。

由图5可知，滚动轴承在525点左右出现早期故障，在525点以前轴承对应的可靠度为1。之后可靠度开始逐渐减小，到860点时滚动轴承可靠度已经降到0.25附近，此时轴承处于故障程度比较严重的状态。

利用IFOA对SVR中的参数C，g以及ε进行寻优，L₀＝30，Sizepop＝20，Maxgen＝100。训练数据取第680点到第860点，输入向量维数为20，建立IFOA-SVR预测模型。表2为利用不同的优化算法对SVR进行寻优得到的参数值。

表2利用不同算法寻优得到的3个参数值

预测数据从第861点开始，共进行了50次迭代运算，通过长期迭代预测法得到最终预测结果。如果只优化SVR模型中的C和g，预测误差会很大并且预测步长较短。

利用GA对SVR模型中的C、g以及ε同时进行寻优，预测结果如图6所示。

利用PSO对SVR模型中的C、g以及ε同时进行寻优，预测结果如图7所示。

利用DS-FOA对SVR模型中的C、g以及ε同时进行寻优，预测结果如图8所示。

利用IFOA对SVR模型中的C、g以及ε同时进行寻优，并将MAE、RMSE、NMSE以及MAPE的和作为适应度函数，预测结果如图9所示。

实际退化曲线在890点以后呈上升趋势，由图6可以看出，预测曲线在890点以后，明显偏离实际退化曲线，预测效果不好；由图7可以看出，预测曲线在890点以后上升趋势不明显；由图8可以看出，预测曲线与实际退化曲线较为接近，但是在893至901点，与实际退化曲线相差较大；对比图6至图9可以看出，IFOA优化SVR所得出的预测曲线与实际性能退化曲线更为接近，预测效果更好。

为了定量地评估预测结果，分别计算所预测50个点的MAE、RMSE、NMSE和MAPE，计算结果如表3所示。

表3 50点预测误差分析

由表3可知，在本发明对SVR中3个参数均优化的情况下，进行50点预测时，与GA、PSO以及DS-FOA优化的SVR预测模型相比，IFOA-SVR模型的4类预测误差均最小，其4类误差和也最小。证明IFOA具有更好的空间寻优能力，能够准确地找到SVR模型中参数的最佳值，体现了IFOA-SVR模型在滚动轴承性能退化趋势预测方面的优势。

与文献[7]到[9]对比分析可知，本发明所得实验结果MAPE低于文献[7]，预测步长增加了30点；与文献[8]相比，MAPE在同一数量级内稍高，但预测步长增加了43点；相比于文献[9]，在降低了MAE和MAPE的同时预测步长增加了20点。

将预测所得50点数据嵌入可靠度模型中，可以计算出这50点的可靠度值，可靠度曲线如图10所示。

从图10中可以看出，第861点的可靠度值在0.25附近，并且这50点的可靠度呈下降趋势，到910点时可靠度已经下降到0.08左右。这表明轴承运行状态很差，已经出现非常严重的故障，应该做好及时更换轴承的准备。

通过对本发明方法的应用得出如下结论：

(1)提出Hilbert包络解调与数学形态学分形维数相结合的特征提取方法。通过控制变量法得出最大分析尺度λ_max以及结构元素Se这两个参数的最佳设置，计算出包络信号的出数学形态学分形维数并作为刻画滚动轴承性能退化过程的特征。

(2)将MAE、RMSE、NMSE以及MAPE的和作为IFOA的适应度函数，对SVR模型中的3个参数进行寻优，建立了IFOA-SVR模型。利用该模型对滚动轴承的性能退化趋势进行了50点预测，与GA，PSO和DS-FOA优化的SVR模型相比预测误差最小。同时，与相关文献进行了横向对比，在保证预测精度的前提下，增加了预测步长。

(3)利用IFOA求解WPHM的待定参数，建立了可靠度模型。将IFOA-SVR模型所得的特征预测结果嵌入到可靠度模型中，得到可靠度预测曲线，能够预测滚动轴承在未来50点的运行状态可靠度。

本发明中涉及的参考文献：

[1]何正嘉,曹宏瑞,訾艳阳,等.机械设备运行可靠性评估的发展与思考[J].机械工程学报,2014,50(2)：171-186.

HE Zhengjia,CAO Hongrui,ZI Yanyang,et al.Developments and Thoughts onOperational Reliability Assessment of Mechanical Equipment[J].Journal ofMechanical Engineering,

2014,50(2)：171-186.

[2]WANG Zhonglai,HUANG Hongzhong,LI Yanfeng.An approach toreliability assessment under degradation and shock process[J].IEEETransactions on Reliability,2011,60(4)：852-863.

[3]何群,李磊,江国乾,等.基于PCA和多变量极限学习机的轴承剩余寿命预测[J].中国机械工程,2014,25(7)：984-988+989.

HE Qun,LI Lei,JIANG Guoqian,et al.Residual Life Predictions forBearings Based on PCA and MELM[J].China Mechanical Engineering,2014,25(7)：984-988+989.

[4]申中杰,陈雪峰,何正嘉,等.基于相对特征和多变量支持向量机的滚动轴承剩余寿命预测[J].机械工程学报,2013,49(2)：183-189.

SHEN Zhongjie,CHEN Xuefeng,HE Zhengjia,et al.Remaining LifePredictions of Rolling Bearing Based on Relative Features and MultivariableSupport Vector Machine[J].Journal of Mechanical Engineering,2013,49(2)：183-189.

[5]陈向民,于德介,李蓉.基于形态分量分析和包络谱的轴承故障诊断[J].中国机械工程,2014,25(8)：1047-1052+1053.

CHEN Xiangmin,YU Dejie,LI Rong.Fault Diagnosis of Rolling BearingsBased on Morphological Componenet Analysis and Envelope Spectrum[J].ChinaMechanical Engineering,2014,25(8)：1047-1052+1053.

[6]张亢,程军圣,杨宇.基于局部均值分解与形态学分形维数的滚动轴承故障诊断方法[J].振动与冲击,2013,32(9)：90-94.

ZHANG Kang,CHENG Junsheng,YANG Yu.Roller bearing fault diagnosisbased on local mean decomposition and morphological fractal dimension[J].Journal of Vibration and Shock,2013,32(9)：90-94.

[7]陈昌,汤宝平,吕中亮.基于威布尔分布及最小二乘支持向量机的滚动轴承退化趋势预测[J].振动与冲击,2014,33(20)：52-56.

CHEN Chang,TANG Baoping,LU Zhongliang.Degradation trend prediction ofrolling bearings based on Weibull distribution and least squares supportvector machine[J].Journal of Vibration and Shock,2014,33(20)：52-56.

[8]宋梅村,蔡琦.基于支持向量回归的设备故障趋势预测[J].原子能科学技术,2011,45(8)：972-976.

SONG Meicun,CAI Qi.Fault Trend Prediction of Device Based on SupportVector Regression[J].Atomic Energy Science and Technology,2011,45(8)：972-976.

[9]朱霄珣.基于支持向量机的旋转机械故障诊断与预测方法研究[D].北京：华北电力大学,2013.

ZHU Xiaoxun.Research on Rotating Machine Fault Diagnosis and StatePrediction Method Based on Support Vector Machine[D].Beijing：North ChinaElectric Power University,2013.

[10]P.J.García Nieto,E.García-Gonzalo,F.Sánchez Lasheras et al.HybridPSO-SVM-based method for forecasting of the remaining useful life foraircraft engines and evaluation of its reliability[J].Reliability Engineeringand System Safety,2015,138：219-231.

[11]PAN W T.A new fruit fly optimization algorithm:Taking thefinancial distress model as an example[J].Knowledge-Based Systems,2012,26(2)：69-74.

[12]宁剑平,王冰,李洪儒,等.递减步长果蝇优化算法及应用[J].深圳大学学报理工版,2014,31(4)：367-373.

NING Jianping,WANG Bing,LI Hongru,et al.Research on and applicationof diminishing step fruit fly optimization algorithm[J].Journal of ShenzhenUniversity Science and Engineering,2014,31(4)：367-373.

[13]张义民,孙志礼.机械产品的可靠性大纲[J].机械工程学报,2014,50(14)：14-20.

ZHANG Yimin,SUN Zhili.The Reliability Syllabus of Mechanical Products[J].Journal of Mechanical Engineering,2014,50(14)：14-20.

[14]丁锋,何正嘉,訾艳阳,等.基于设备状态振动特征的比例故障率模型可靠性评估[J].机械工程学报,2009,45(12)：89-94.

DING Feng,HE Zhengjia,ZI Yanyang,et al.Reliability Assessment Basedon Equipment Condition Vibration Feature Using Proportional Hazards Model[J].Journal of Mechanical Engineering,2009,45(12)：89-94.

[15]陈昌.基于状态振动特征的空间滚动轴承可靠性评估方法研究[D].重庆：重庆大学,2014.

CHEN Chang.Reliability Assessment Method for Space Rolling BearingBased on Condition Vibration Feature[D].Chongqing：Chongqing University,2014.

[16]董绍江.基于优化支持向量机的空间滚动轴承寿命预测方法研究[D].重庆：重庆大学,2012.

DONG Shaojiang.Research on Space Bearing Life Prediction Method Basedon Optimized Support Vector Machine[D].Chongqing：Chongqing University,2012.

[17]Rexnord technical services,bearing data set,IMS,University ofCincinnati,NASA Ames Prognostics Data Repository,NASA Ames,Moffett Field,CA,2007[2014-04].http://ti.arc.nasa.gov/project/prognostics-data-repository.

Claims (2)

1.一种基于数学形态学和IFOA-SVR的滚动轴承可靠度预测方法，其特征在于：所述方法的实现过程为：

步骤一、求出包络信号的数学形态学分形维数D_M；

步骤二、从获取的数学形态学分形维数D_M选取训练样本对，基于所述训练样本对利用IFOA对SVR模型中的参数C、g以及ε同时进行寻优，建立IFOA-SVR预测模型：

构建IFOA的过程：

(1)初始化算法参数：设置果蝇种群规模Sizepop，最大觅食代数Maxgen，并随机初始化果蝇群体位置坐标(X₀,Y₀,Z₀)；

(2)果蝇个体利用嗅觉随机搜索的方向和距离可以通过式(12)获得

$\left\{\begin{array}{c}{X}_i=X_0+(L_0-\frac{L_0(gen-1)}{Maxgen})\&times;rands(1,1)\\ Y_i=Y_0+(L_0-\frac{L_0(gen-1)}{Maxgen})\&times;rands(1,1)\\ Z_i=Z_0+(L_0-\frac{L_0(gen-1)}{Maxgen})\&times;rands(1,1)\end{array}\right.---\left(12\right)$

式中，i＝1,2,...,Sizepop，L₀为初始步长值，gen为当前觅食代数；

(3)由于无法确定食物源的具体位置，所以需要通过式(13)估计第i个果蝇个体的当前位置与坐标原点间的距离Dist_i，之后计算出味道浓度判定值S_i

${\mathrm{Dist}}_i=\sqrt{{X_i}^2+{Y_i}^2+Z_i^2}---\left(13\right)$

S_i＝1/Dist_i (14)

(4)将S_i代入味道浓度判定函数，计算出果蝇个体当前位置的味道浓度

Smell_i＝function(S_i) (15)

(5)当前果蝇群体中具有最高味道浓度的个体，可由式(16)获得

[bestSmell,bestIndex]＝max(Smell_i) (16)

式中，bestSmell表示果蝇群体中具有最高味道浓度的个体的味道浓度值，bestIndex表示果蝇群体中具有最高味道浓度的个体的位置；

(6)保留果蝇群体中最佳味道浓度值和与其对应的个体坐标，同时果蝇群体利用自身的视觉对食物源进行定位，然后飞往食物源所在的位置；

$\{\begin{array}{c}Smellbest=bestSmell\\ X_0=X\left(bestIndex\right)\\ Y_0=Y\left(bestIndex\right)\\ Z_0=Z\left(bestIndex\right)\end{array}---\left(17\right)$

(7)进入迭代寻优过程，重复步骤(2)-(5)，并判断当前味最高道浓度是否好于前一迭代味道浓度，且gen<Maxgen；若成立，则执行步骤(6)；

建立IFOA-SVR模型，其过程为：

(1)初始化IFOA参数，包括Sizepop、Maxgen以及果蝇个体初始位置；对SVR中的3个参数：惩罚系数C、核函数参数g以及不敏感误差ε进行寻优，所以初始坐标为(X₀ ¹,Y₀ ¹,Z₀ ¹)，(X₀ ²,Y₀ ²,Z₀ ²)和(X₀ ³,Y₀ ³,Z₀ ³)；

(2)附与每个果蝇个体随机方向和飞行距离，并用递减搜索步长代替固定步长，得到(X_i ¹,Y_i ¹,Z_i ¹)，(X_i ²,Y_i ²,Z_i ²)，(X_i ³,Y_i ³,Z_i ³)，计算当前果蝇个体与原点间距离的倒数，得到味道浓度判定值S_i ¹，S_i ²和S_i ³；

(3)确定SVR中参数C、g和ε的范围，即C∈[2^-14,2¹⁴]，g∈[2^-14,2¹⁴]，ε∈[2^-14,2¹⁴]；

(4)将所述训练样本对输入到SVR模型中，对模型进行训练，将平均绝对误差MAE、均方根误差RMSE、归一化均方误差NMSE以及平均绝对百分误差MAPE的和作为适应度函数，即

$\begin{array}{c}{\mathrm{Smell}}_i=Fitness(C_i,g_i,{\mathrm{\&epsiv;}}_i)\\ =\&lsqb;MAE+RMSE+NMSE+MAPE\&rsqb;(C_i,g_i,{\mathrm{\&epsiv;}}_i)\end{array}---\left(18\right)$

(5)找到适应度函数Fitness的最小值对应的果蝇个体，开始迭代寻优，并判断最小Fitness是否低于前一代最小Fitness；如低于，则保留最小Fitness值及其对应的坐标，并将其赋给初始坐标；如高于，则返回步骤(2)；

(6)找到C、g和ε的最佳值，建立IFOA-SVR预测模型；

步骤三、从所述性能退化状态特征中截取一定长度Z作为WPHM的响应协变量，结合极大似然估计获得似然函数方程组，把方程组中每个方程绝对值的和作为IFOA的适应度函数，求出方程组的解，确定WPHM的待定参数，进而得到可靠度模型R(t,Z)，过程如下：

基于WPHM建立滚动轴承性能退化状态特征与可靠度之间的数学关系，WPHM的表达式为

$h(t,Z)=\left(\frac{\mathrm{\&beta;}}{\mathrm{\&eta;}}\right){\left(\frac{t}{\mathrm{\&eta;}}\right)}^{\mathrm{\&beta;}-1}\exp \left(\mathrm{\&mu;}Z\right)---\left(19\right)$

式中β为形状参数，η为尺度参数，μ为协变量回归参数，t表示时间，Z为截取一定长度的性能退化状态特征，即从数学形态学分形维数D_M中截取一定长度数据；

h(t,Z)和可靠度函数R(t,Z)之间的关系为

$h(t,Z)=-\frac{d}{dt}lnR(t,Z)---\left(20\right)$

从而，R(t,Z)可表示为

$\begin{array}{c}R(t,Z)=\exp (-{\&Integral;}_0^{t}h(t,Z\left)\right)\\ =\exp (-{\&Integral;}_0^t\frac{\mathrm{\&beta;}}{\mathrm{\&eta;}}{\left(\frac{t}{\mathrm{\&eta;}}\right)}^{\mathrm{\&beta;}-1}\exp (\mathrm{\&mu;}Z\left)dt\right)\\ =\exp (-{\left(\frac{t}{\mathrm{\&eta;}}\right)}^{\mathrm{\&beta;}}\exp (\mathrm{\&mu;}Z\left)\right)\end{array}---\left(21\right)$

利用极大似然估计法得到关于待定参数β、η、μ的方程组并利用IFOA对该方程组进行求解，确定WPHM的表达式h(t,Z)的待定参数β、η、μ的最终取值；将参数β、η、μ代入式(21)即可建立可靠度模型；

步骤四、将性能退化状态特征的最后一部分数据作为IFOA-SVR预测模型的输入，采用长期迭代预测法获取性能退化状态特征预测结果；

步骤五、将IFOA-SVR预测模型得出的预测结果作为响应协变量嵌入到可靠度模型R(t,Z)中，即可计算出所述性能退化状态特征预测结果所对应的可靠度，实现滚动轴承可靠度的预测。

2.根据权利要求1所述一种基于数学形态学和IFOA-SVR的滚动轴承可靠度预测方法，其特征在于：

在步骤一中，通过Hilbert变换计算出状态监测数据的包络信号；

步骤一一、获取滚动轴承的状态监测数据：滚动轴承的振动信号；

步骤一二、通过Hilbert变换计算出状态监测数据的包络信号；

步骤一三、利用式(4)求出每次采集振动信号的包络信号的数学形态学分形维数D_M，将D_M作为滚动轴承性能退化状态特征；

$log\left(\frac{A_{Se}\left(\mathrm{\&lambda;}\right)}{{\mathrm{\&lambda;}}^2}\right)=D_{M}log\left(\frac{1}{\mathrm{\&lambda;}}\right)+c,\mathrm{\&lambda;}=1,2,...,{\mathrm{\&lambda;}}_{max}---\left(4\right)$

式中：D_M为信号的Minkowski-Bouligand维数即数学形态学分形维数，c为常数，对和log(1/λ)进行最小二乘线性拟合即可得到对信号数学形态学分形维数D_M；λ＝1,2,…,λ_max，λ_max为最大分析尺度；

A_Se(λ)为定义在分析尺度λ下f(n)关于Se(n)膨胀和腐蚀运算的覆盖面积，f(n)表示滚动轴承振动信号的包络信号；Se(n)为结构元素，为一维离散向量；Se(n)的长度及λ_max的值通过控制变量法求得。