轴承性能时间序列的概率信息求取及退化分析*

常 振，魏剑波，平晓明，曹茂来，王二化

(1.杭州轴承试验研究中心有限公司，浙江 杭州 310022；2.机械工业轴承产品质量检测中心(杭州)， 浙江 杭州 310022；3.常州信息技术学院 智能装备学院，江苏 常州 213164)

0 引 言

轴承性能时间序列是将某种性能属性的数值按照时间测量的先后顺序排列所形成的数列，根据时间序列所反映出来的发展过程、方向和趋势，进行类推或延伸，借以预测或分析下一段时间或若干年后可达到的动态运转水平，进而快速有效地进行产品性能故障诊断[1-3]。

时间序列的信息挖掘一直受到学术界与工程界的高度重视，特别是在航空航天、金融、经济、天文、地质等应用十分广泛[4-7]。而早期时间序列分析法的模型几乎都是线性的，随产品质量指标的要求不断提高，研究人员发现非线性模型才能合理解释工程实际问题[8，9]。轴承服役过程中伴随有振动、温度、摩擦力矩等不同属性的非线性性能时间序列，因此，可从该类时间序列中提取系统参与演化的所有变量的大量信息，从而用于轴承某些性能信号概率信息的求取与退化特征分析。

轴承性能时间序列概率密度函数的求取是数据处理与信号分析的根基，可据此获取某一性能信号的真值估计与区间波动。

本文基于轴承的非线性性能时间序列，首先通过自助法建立时间序列仿真数据，运用最大熵原理对大量的仿真数据进行概率密度求取，分析其区间波动情况并验证模型的准确性；然后基于模糊集合理论[10，11]进行轴承性能时间序列的退化特征分析；两套模型兼并使用综合讨论轴承性能信号的演变规律并识别其性能退化历程。

1 轴承性能时间序列的概率信息求取

1.1 性能时间序列分组处理

设轴承某一性能时间序列向量表示为：

X=(x(1),x(2),…,x(n),…x(N))

(1)

式中：x(n)—轴承性能时间序列中的第n个数据，n=1,2,…,N；N—样本个数。

取前n个数据序列作为训练组，记为XA；后N-n个数据序列作为验证组，记为XB，即该轴承的性能时间序列分组为训练组XA和和验证组XB，表示为：

XA=(x(1),x(2),…,x(n))

(2)

XB=(x(n+1),x(n+2),…,x(N))

(3)

1.2 训练组XA灰自助再抽样

为满足灰预报模型GM(1,1)关于x(n)≥0的苛刻要求，在式(2)中，若有x(n)<0，则人为选取一个常数c，使得x(n)+c≥0即可。

实际分析时，X表示为：

XA=(x(n)+c;t=1,2,…,N)

(4)

从XA中取与时刻n紧邻的前m个时刻的数据(包括时刻n的数据)，构成时刻n的动态分析子向量，即：

Xm=[xm(u)];u=n-m+1,n-m+2,…,n;n≥m

(5)

运用自助法，在时刻n，从Xm中等概率可放回地随机抽取一个数，共抽取m次，可得到一个自助样本Y1，它有m个数据。

按此方法重复执行B次，得到B个样本，可表示为：

YBootstrap=(Y1,Y2,…,Yb,…,YB)

(6)

式中：Yb—第b个自助样本；B—总的自助再抽样次数，也是自助样本的个数。

且有：

Yb=[yb(u)]

(7)

其中：u=n-m+1,n-m+2,…,n；b=1,2,…,B。

根据灰预报模型GM(1,1)[12]，可得到ω=n+1时刻，组成样本容量为B的大数据模拟信号XGBootstrap，可表示为如下向量：

(8)

1.3 基于最大熵原理求解模拟信号XGBootstrap的概率密度函数

(9)

随后可根据最大熵方法获得基于模拟信号XGBootstrap的密度函数的最优估计。最大熵方法的核心原理是在所有的可行解中，满足熵最大的解是最“无偏”的，亦即：

(10)

在具体的求解过程中，引入拉格朗日乘子，可以把概率分布求解问题转化为对拉格朗日乘子的求解问题[13]。

1.4 积分区间的映射求解概率密度函数

为了使计算过程的收敛，笔者将序列XGBootstrap中的数据按进行递增排序并分成β组，随后画出对应的直方图。同时，可得到组中值zμ和频数Γμ，u=2,3,…,β+1。然后将直方图扩展成β+2组，并令Γ1=Γξ+2。将原始数据区间S映射到区间[-e,e]中，e—自然常数。

令：

(11)

(12)

(13)

b=e-azξ+2

(14)

式中：e—自然常数，e=2.718 282。

因此，所求的最大熵原理构建的轴承性能时间序列概率密度函数变换为：

(15)

式中：j—原点矩阶数，常用j=5；c0—首个拉格朗日乘子；ci—第i+1个拉格朗日乘子，i=1,2,…,j。

对于实数α双侧分位数，有概率：

(16)

(17)

式中：XU，XL—轴承性能时间序列估计区间的上界值和下界值；[XL,XU]—估计区间。

1.5 轴承性能时间序列概率信息求取方法的验证

已知该轴承性能时间序列的前n个数据序列作为训练组XA，后N-n个数据序列作为验证组XB。其中利用训练组XA获取上述性能时间序列的概率密度函数和对应的区间估计，则用验证组XB来证实上述方案的可行性：

假设在N-n个验证组XB的数据中有y个数据在估计区间[XL,XU]之内，则上述模型的可行度或准确度表示为：

(18)

此处所提准确度也是指验证组XB中的数据落入估计区间的频率，实现了所提模型的自我验证；同时也获取了该轴承性能时间序列的区间上、下限，可以为工程中及时掌控轴承的服役工况提供理论支持。

2 轴承性能时间序列的退化特征分析

2.1 性能时间序列向模糊关系的转换

轴承性能退化特征分析时，使用的是模糊等价关系转换。然而，在工程实际中往往得到的是模糊相似关系，必须将这种模糊相似关系转变为模糊等价关系。传递闭包法可以进行这种关系的转换，以获得轴承性能时间序列的模糊等价关系，具体求解步骤与方法如下：

将轴承性能时间序列分组处理：平均分为p组，设每组有q个数据，即一小组可组成数据序列：

Zl=(Zl1,Zl2,…,Zlk,…,Zlq)k=1,2,…,q,l=1,2,…,p

(19)

式中：p—样本数；q—每个样本的样本量；Zlk—第l个样本的第k个数据。

则该轴承性能时间序列构成一个新的集合，即：

(20)

对于任意的模糊关系R，如若存在如下关系：

T(R)=Rh-1=Rh=…h=1,2,3,…

(21)

式中：T(R)—模糊数学中的模糊等价关系。

T(R)可以按照上述式子逐步求出，具体过程如下：

第一步，求出R2=R∘R；

第二步，求出R4=R2∘R2；

……

第u步，直到求出R2u=Ru为止。

其中：运算符“°”—矩阵的模糊运算M(∧,∨)；“∨”—“或”运算取最大；“∧”—“与”运算取最小；Ru—所求的该轴承性能时间序列模糊等价关系T(R)，也是模糊集合理论的传递闭包。

即有：

T(R)=Ru

(22)

有：

(23)

其中：0≤αil≤1，且满足：

(24)

式中：αil—性能时间序列分组后第i个样本Zi和第l个样本Zl两者之间的模糊等价关系，亦即Zi样本信息和Zl样本信息的符合程度。

2.2 轴承性能退化依据

在工学、文学、哲学等各大领域中，0和1都可以表示事物的真假、有无、强弱等两个极端状态。然而任何事物都有着或多或少的关联，很难用0和1的确定性规则判定两个事物的必然规律。

可以在模糊等价关系αil中设定λ阈值来诊断两者关系密切的大小或者产生特征退变存在的显著性：

(1)若αil>λ，则两段轴承的性能时间序列Zi和Zl在λ阈值水平下彼此之间的相似度很高，关联程度较大，即表明两个系统间不存在性能退化；

(2)若αil≤λ，则两段轴承的性能时间序列Zi和Zl在λ阈值水平下彼此之间的相似度很低，关联程度较小，即表明两个系统间变化大、演变强烈，存在有性能退化趋势。

根据模糊集合理论，λ为研究对象从一个极端属性过渡到另一极端属性的边界，也是判定两段轴承性能时间序列是否有退变产生的显著性。在工程实际中，当λ=0.5时两系统之间的模糊关系达到最大，介于最难分辨的真假、有无、强弱之间；当λ>0.5时两系统间关系趋于清晰，相似度较高且可判定两者未发生退变[14，15]。

所以在应用过程中，取λ=0.5来界定轴承性能是否发生特征退化。

2.3 轴承性能退化的模糊表征

定义退化系数集合为：

U=(u1,u2,…,uξ,…,uw-1)

(25)

其中，有：

(26)

式中：uξ—轴承性能系数，即模糊等价关系αil的分段均值；w—样本含量；ξ—各样本采样的时间先后顺序，是一个时间变量。

根据退化系数uξ来判断轴承性能时间序列的演变历程：(1)如若uξ>λ=0.5，则表明轴承服役期间保持良好的运转特性，性能未发生而恶性退化；(2)uξ<λ=0.5，则表明其服役期间各段性能数据相似度下降，退变趋势明显，性能发生明显退化甚至轴承已发生损坏。

3 实验及结果分析

案例一：滚动轴承振动性能时间序列试验研究。

该试验在ABLT-1型轴承寿命强化试验机上进行的，试验轴承为6 208。转速为4 000 r/min，纯径向载荷Fr=8 kN，油润滑。每10 min采集振动均方根值一次，单位ms-2，共采集60个数据，获得振动性能时间数据序列X振动，如图1所示。

图1 轴承振动性能时间数据序列X振动

案例二：滚动轴承温度性能时间序列试验研究。

该试验在ABLT-3型轴承寿命强化试验机上进行的，试验轴承为6 203。转速为7 000 r/min，纯径向载荷Fr=2.5 kN，油润滑。每10 min采集温度值一次，单位℃，共采集60个数据，获得温度性能时间数据序列X温度，如图2所示。

图2 轴承温度性能时间数据序列X温度

案例三：滚动轴承摩擦力矩性能时间序列试验研究。

研究对象为轴承组件的稳态电流信号(单位：mA)，反作用控制箱输出指令电压带动真空实验装置中的滚动轴承转动，由反馈装置取样并转换后将得到的电流信号反馈给控制箱。反馈得到的稳态电流信号间接得到滚动轴承运转期间摩擦力矩时间序列。等间隔地采集轴承摩擦力矩60次，获得摩擦力矩性能时间序列X摩擦力矩，如图3所示。

图3 轴承摩擦力矩性能时间数据序列X摩擦力矩

3.1 概率密度求取

案例中性能时间序列X，笔者取前10个数据作为训练组，记为XA。从训练组XA序列中等概率可放回地进行抽样，每次抽样个数v=6；连续上述步骤重复抽取B=10 000次，进而可组成样本容量为10 000的大数据模拟信号XGBootstrap。

案例一生成的灰自助数据如图4所示。

图4 振动性能信号生成的灰自助数据XGBootstrap

案例二生成的灰自助数据如图5所示。

图5 温度性能信号生成的灰自助数据XGBootstrap

案例三生成的灰自助数据如图6所示。

图6 摩擦力矩性能信号生成的灰自助数据XGBootstrap

案例一求得的概率密度函数如图7所示。

图7 振动性能信号概率密度函数

案例二求得的概率密度函数如图8所示。

图8 温度性能信号概率密度函数

案例三求得的概率密度函数如图9所示。

图9 摩擦力矩性能信号概率密度函数

3.2 区间预报

根据图(7～9)概率密度函数信息(由前n=10个数据XA获得)，并给定α=0.002 7显著性水平下可获得轴承性能信号的区间估计，相应置信水平为99.73%。用剩余的60-n=50个数据来检验该预报的可靠性，并掌控后续信号的退变信息。

其结果如表1所示。

表1 时间序列X振动、X温度、X摩擦力矩的区间预报分析

据上述滚动轴承性能时间序列概率密度信息的求取过程可知，数据处理过程中仅用了训练组的10个数据，且对该组数据无任何先验信息也无任何假设约束条件，表明该模型具有优越的客观主动性；

随后以验证组来证明其区间预报的准确性，在99.73%的置信水平下，使用灰自助最大熵法区间估计的上、下边界差值较小，其中，性能时间序列X振动的区间差值为0.42；性能时间序列X温度的区间差值为2.34；性能时间序列X摩擦力矩的区间差值为23.42，同时表明了所提方法具有较高的预报精度。更何况滚动轴承振动、温度、摩擦力矩性能时间序列的随机性很强，且具有明显的不确定度，其区间预报十分复杂，所以该方法的区间预测已相当可观。

此外，在验证组的50个数据进行方法验证过程中，试验数据X振动、X温度中验证组落在预报区间之外的个数极少，分别为2个和1个，表明预报结果的准确可行性，可用于工程实践的在线监测与诊断。然而，试验数据X摩擦力矩中验证组落在预报区间之外的个数很多，为33个，并非表明预报结果不准确可行，而是表明验证组的后50个数据相对于训练组的前10个数据有着明显的差异性，亦即轴承服役过程中随时间的演变，其摩擦力矩性能信号产生了严重的退化。由图3的原始数据X摩擦力矩可知，摩擦力矩信号呈现上升趋势，即表明轴承运行状况并非稳定，内部摩擦逐渐增大，其性能信号已发生变异或退化。

基于灰自助最大熵法X振动、X温度的区间误报率分别为4%和2%，即表明此次性能时间序列预报的可靠度在96%～98%，所以该方法进行时间序列区间的预报精度高且可靠性强；同时也说明后50个验证组组数据与前10训练组数据保持着良好的一致性，退化迹象十分不明显，即轴承运转期间保持较高的稳定性和平稳性。

基于灰自助最大熵法X摩擦力矩的区间误报率66%，即表明此次性能时间序列预报的可靠度在44%，结果表明后50个验证组组数据与前10训练组数据一致性不强，退化迹象较为明显，即轴承运转期间已产生明显退化和变异。

3.3 轴承退化性能分析

案例一。

数据处理后得到的模糊等价关系矩阵T(R)为：

(27)

稳定系数集合U=[0.71 0.72 0.70 0.72 0.70 0.71 0.68 0.66 0.60]，结果如图10所示。

案例二。

数据处理后得到的模糊等价关系矩阵T(R)为：

图10 轴承振动性能时间序列X振动退化系数

(28)

稳定系数集合U=[0.64 0.61 0.61 0.61 0.63 0.61 0.61 0.61 0.61]，结果如图11所示。

案例三。

数据处理后得到的模糊等价关系矩阵T(R)为：

图11 轴承温度性能时间序列X温度退化系数

(29)

稳定系数集合U=[0.70 0.67 0.64 0.62 0.59 0.56 0.51 0.50 0.48]，结果如图12所示。

图12 轴承摩擦力矩性能时间序列X摩擦力矩退化系数

由上述轴承性能时间序列X振动、X温度、X摩擦力矩的3个案例的退化信息处理结果可知：

案例一、二的振动性能和温度性能表现稳定，最小退化系数分别为uξ=0.600和uξ=0.609，两者均大于λ=0.5阈值；则表明案例一的轴承振动性能和案例二的温度性能未发生明显异常，前后时间序列未发生明显变异，轴承退化特征不明显，轴承运转状况稳定可靠，健康状况较为良好；

案例三的摩擦力矩性能表现不良，最小退化系数为uξ=0.477<0.5阈值，则表明案例三轴承的摩擦力矩信号已发生明显异常，前后时间序列有明显变异，轴承性能信号已发生明显退化迹象，图12中退化系数呈现出明显的减小趋势，该结果与其原始数据系列逐渐增大有关，且与3.2部分区间预报分析结果保持良好的一致性，均说明了该轴承摩擦力矩性能信号退化明显，运转状况不再安全可靠；为避免恶性事故的发生，应及时维护或更换轴承，并对其进一步健康监控。

综上可知，轴承性能时间序列的概率信息的求取为更深层次挖掘其隐含信息奠定了基础，据此进行的区间预测实现了模型准确性的自我验证，并可以较好地反映其工作期间性能信号的波动情况；另外，轴承性能时间序列的退化系数，可以很好地识别出其特定性能的演变规律和退化趋势。

因此，笔者所提出的自助最大熵结合法可以有效地拟合轴承性能时间序列的概率密度函数，并可预测其性能区间波动状态；基于模糊理论的退化系数可准确地监控其性能退化动态信息。

4 结束语

针对轴承性能时间序列概率信息求取及退化分析问题，笔者对轴承振动、温度、摩擦力矩3种性能时间序列进行了研究，结论如下：

(1)在概率信息求取过程中，X振动、X温度区间误报率低，概率信息求取较为准确可靠；X摩擦力矩的误报率为66%，这主要是因为该性能信号具有明显的退化信息，而非所提模型不可靠；

(2)在退化性能分析过程中，X振动、X温度表现出良好的服役状况；X摩擦力矩的退化系数uξ=0.477<0.5，具有明显退化迹象；与概率信息求取的结果保持较好的一致性；

(3)结合轴承性能时间序列的概率信息和退化系数，可有效评估轴承性能时间序列的变化趋势本质特征和退化演变迹象，实时掌控轴承的健康状况，并及时发现其失效隐患。