卜德天，陈龙，刘红彬，吕文静，张弘毅
(1.河南科技大学 机电工程学院，河南 洛阳 471003；
2.洛阳LYC汽车轴承科技有限公司，河南 洛阳 471003；
3.洛阳轴承研究所有限公司, 河南 洛阳 471039)
汽车轮毂轴承[1](Hub Bearing)是轮毂总成的重要组成部分，其疲劳寿命与可靠性计算的确定对汽车行驶安全性有重大影响。Lundberg和Palmgren基于赫兹接触理论和威布尔材料强度统计理论建立了滚动轴承疲劳寿命估计方法，被广泛接受并成为轴承行业的计算标准[2]。
Gustafsson[3]于1962年首次提出采用加速度传感器采集轴承振动信号分析的方法，振动监测随后广泛应用于回转机械的运行状态监测。相比于基于声信号[4]、静电信号[5]和定子电流[6]的监测方法，振动监测不仅能够形象、直观地描述轴承性能，还能够从振动信号中提取性能退化过程的特征信息[7-8]。振动监测获取的大量状态演变数据也为基于振动的轴承性能演变以及剩余寿命预测提供了充足的数据源。
依据振动数据测量时间的先后关系，一些学者采用振动时间序列方法开展机械设备的运行状态预测：文献[9]以神经网络为基础构建振动频率时间序列预测方法，提高了旋转机械设备运行状态预测精度；
文献[10]开展潮流能水平轴水轮机振动位移时间序列的自相似性研究，为采用自相似理论研究潮流能水轮机的动态特性奠定了数学基础，所揭示的时间序列整体与局部之间的关系及其自相似性为潮流能水轮机的故障诊断提供了技术支撑；
文献[11]进行变压器油箱表面振动信号时间序列的混沌特性分析，证明变压器表面振动信号的时间序列具有混沌特性，为变压器振动信号进一步的处理及应用提供了参考；
文献[12]基于时间序列分析机车振动信号，建模和预测结果表明时间序列分析技术能够快速有效地处理振动信号的预测问题。
由于数据量充足，统计分析方法也被广泛应用于分析滚动轴承的振动数据：文献[13]在研究三参数威布尔分布参数的分位数估计、拟估计和改良极大似然估计时采用随机模拟法对其优良性进行了探究，并给出了可靠度的置信区间；
文献[14]在研究叶片发动机的故障实例时将最小二乘法和平均秩次法结合对威布尔参数进行计算；
文献[15]基于Nelson定律，采用等累积退化量和等累积失效概率的模型得到了不同应力下的等效试验时间；
文献[16]在研究铁路设备关键零部件的可靠性分析中提出了一种新的相关系数法，对收集数据进行三参数威布尔分布参数估计并建立了铁路设备关键零部件可靠性分析模型；
文献[17]在对导弹武器系统的性能检测时研究了小样本情况下威布尔参数的计算；
文献[18]以铰销式碳纤维增强基复合材料环带为研究对象，基于SWT疲劳寿命预测模型推导了服从威布尔分布的概率疲劳寿命预测方法并证明了此预测模型的正确性。
本文在研究汽车轮毂轴承疲劳寿命试验振动数据时，将振动时间序列划分为多个子序列，对各子序列的数据开展统计分析，以研究不同序列中威布尔分布参数的变化，总结形状参数与尺度参数的变化规律，实现轴承疲劳失效的提前预报。
试验轴承为双列角接触汽车轮毂轴承，型号为DAC2F11410561，内径为34 mm、外径为83 mm，球材料为GCr15，球直径为12.7 mm，球组节圆直径为57 mm，接触角为35°。
1.1 试验设备
采用SYJ-LG-NJ轮毂轴承耐久试验机进行试验，如图1所示，试验机由试验主体部件、床身部件、传动系统、加载系统、润滑系统、电气控制系统、计算机软件系统等组成。试验机采用计算机自动控制，可按预设试验载荷、转速和运行时间等主要参数自动运行，在试验过程中可实时获取轴承每个时间点的转速、电流、振动加速度、载荷、油温、轴温等数据，试验机参数见表1。
1—轴向油缸；
2—径向油缸；
3—试验轴承；
4—主轴；
5—电动机。图1 试验机结构示意图Fig.1 Structure diagram of tester
表1 试验机参数Tab.1 Parameters of tester
1.2 试验工况
试验为定速、定载加速寿命试验，径向载荷为8.8 kN，轴向载荷为5.4 kN，转速为400 r/min。试验结束后轴承单元状态和轴承内圈沟道状态如图2所示。
(a) 轴承单元 (b) 轴承内圈沟道图2 试验结束后轴承单元和轴承内圈沟道的状态Fig.2 Status of bearing element and bearing inner ring raceway after test
1.3 试验数据
按照现有疲劳寿命试验经验，试验过程中被测轴承振动加速度值达到初始值2倍时判定轴承发生疲劳失效，试验终止。本试验历时8 320 min，试验过程中每5 s提取一个振动加速度值，总计提取99 840个数据。
时间序列是指某一变量或指标的数值或观测值按时间顺序以相同时间间隔排列的一组数据。为表征轴承振动的演变状态，将振动数据按测量时间先后等间隔排列，构成的振动时间序列如图3所示。
图3 DAC2F1141056轴承振动时间序列Fig.3 Vibration time series of DAC2F1141056 bearing
振动时间序列包含的大量数据为基于统计理论的分析研究提供了基础条件，但一般统计分析方法难以描述时间先后特征。为表征被测轴承试验中的状态演变，将整体时间序列划分为具有时间先后关系的子序列，分析各子序列威布尔分布的形状参数与尺度参数的变化状态，与各子序列对应的振动时间序列比较，实现被测轴承疲劳寿命的威布尔估计。
2.1 分析方法
假设所记录数据(轴承振动值)的总量为N，则轴承振动时间序列Q可描述为
Q=[q(1),q(2),…,q(s),…,q(N)];
s=1,2,…,N，
(1)
式中：q为试验记录的振动值；
s为数据序号。
所有试验数据构成总体样本，依据总体样本的时间先后关系将振动时间序列分为J组，每个子序列包含的样本数量为N/J，子序列可表示为
Qj=[q(j+1),q(j+2),…,q(j+s),…,
q(j+N/J)]。
(2)
不同区间组的统计参数表征被测轴承在统计区间的性能状态，对于第j个子序列，威布尔分布的概率密度函数为
(3)
式中：k，c分别为威布尔分布的形状参数和尺度参数。
第j个子序列的威布尔分布的累积密度函数为
(4)
可改写为
(5)
对等式两边同时取2次对数可得
(6)
y=kx-klgc。
(7)
以上分析过程中可获取不同子序列的形状参数与尺度参数，将不同子序列的统计参数按时间先后关系排序，则可通过参数变化观测轴承状态演变过程的特征。
2.2 分析结果
将99 840个数据分为20组，即J=20，如图3所示。前19组中每组包含5 000个振动数据样本，第20组为剩余的4 840个数据。分别计算20组数据的形状参数和尺度参数，结果见表2：前19个子序列的形状参数和尺度参数在一定范围内波动，无显著规律，说明前19个子序列内轴承状态无明显变化；
第20个子序列内轴承状态则发生了显著变化。
表2 不同子序列的形状参数和尺度参数Tab.2 Shape parameters and scale parameters of different sub-series
依据(3)式、(4)式、(7)式分别绘制各子序列的概率密度函数图、累积分布函数图、直线方程图，结果如图4所示：前19个子序列对应的概率密度函数与累积分布函数分布均较为集中，无显著变化；
第20个子序列对应函数的横坐标跨度明显增加，峰值右移，表征第20个子序列内轴承振动加速度增加，与试验中振动时间序列吻合；
前19个子序列的斜率与截距也无明显变化，而第20个子序列的斜率显著降低，截距显著增加。
以上多个方面的结果均表明：威布尔分布的形状参数和尺度参数能描述轴承状态演变过程，不同区间数据的统计规律能表征轴承状态演变的内在属性。
(a) 概率密度函数 (b) 累积分布函数 (c) 直线方程图4 不同子序列的威布尔分布函数Fig.4 Weibull distribution functions of different sub-series
以上研究采用的是整个疲劳寿命试验的完整寿命周期数据，包含了最终失效数据，而实际工程应用中需要提前预报轴承的疲劳失效，为提升本研究的工程应用价值，选取发生疲劳失效前一段时间内的数据进一步分析。
依据原始试验数据，第98 407个振动加速度为初始振动加速度的2倍，为达到预测目的，截取第88 288～98 287区间内的10 000个数据开展分析，选取的数据区间为由失效值往前120个数据(换算时间单位为10 min)作为研究区间的最后一个数据，截取数据在整个振动时间序列的位置标注(图3)。将所截取数据区间内的10 000个数据划分为10个子序列，每个子序列包含1 000个数据样本。依据上述处理方法得到截取时间段内不同子序列的形状参数和尺度参数，结果见表3，各子序列的威布尔分布函数如图5所示。
表3 截取时间段内的不同子序列的形状参数和尺度参数Tab.3 Shape parameters and scale parameters of different sub-series in interception time period
(a) 概率密度函数 (b) 累积分布函数 (c) 直线方程图5 截取时间段内不同子序列的威布尔分布函数Fig.5 Weibull distribution functions of different sub-series in interception time period
分析可知：前8个子序列形状参数的平均值为24.840 0且数值均在25左右波动，最大值、最小值分别为27.324 7和21.094 5，与均值差异的幅度分别为9.9%和15.1%；
前8个子序列尺度参数的平均值为0.047 9且数值均在0.048左右波动，最大值、最小值分别为0.048 2和0.047 3，与均值差异的幅度分别为0.6%和1.3%；
第9个子序列的形状参数迅速衰减到15.380 6，与前8个子序列均值的变化幅度为38.1%，尺度参数也迅速变化为0.050 2，与前8个子序列均值的变化幅度为4.8%；
第10个子序列的形状参数进一步衰减到8.387 4，与前8个子序列均值的变化幅度达66.3%，尺度参数也进一步变化为0.055 9，与前8个子序列均值的变化幅度为16.7%；
图5中不同子序列威布尔分布函数的对比也清晰表明了第9，10个子序列的显著变化。
第10个子序列最后1个数据换算为时间单位是失效前10 min，假定第10个子序列形状参数或尺度参数的变化幅度达到判定标准，即提前10 min实现被测轴承疲劳失效的预报；
第10个子序列包含1 000个数据，换算为时间单位是83.3 min，假定第9个子序列形状参数或尺度参数变化幅度达到判定标准，则提前93.3 min实现被测轴承疲劳失效的预报。
为验证以上分析方法，选取相同型号、相同试验条件下的另1套DAC2F11410561轴承和相同试验机上试验的另1套不同型号的汽车轮毂轴承的试验数据开展对比验证。
4.1 同型号轴承的验证
另1套DAC2F11410561轴承的疲劳寿命试验记录了136 260个振动加速度值，其振动时间序列如图6所示。采用前述分析方法，以振动加速度达到初始振动加速度的2倍为节点，提前10 min作为分析区间的最后1个数据，选取该轴承100 265～110 264区间内的10 000个数据依次分为10组进行分析，用于验证的数据在整个时间序列的位置如图6所示。
10个子序列的形状参数和尺度参数的变化趋势如图7所示：2种参数的变化规律与表3具有较好的一致性，即前8个子序列的形状参数和尺度参数均在一定范围内波动，第9，10个子序列则发生显著波动，变化趋势和规律与前述分析结果一致。
图6 验证用DAC2F1141056轴承振动时间序列Fig.6 Vibration time series of DAC2F1141056 bearing used for verification
(a) 形状参数 (b) 尺度参数图7 验证用DAC2F1141056轴承截取时间段内不同子序列的形状参数和尺度参数Fig.7 Shape parameters and scale parameters of different sub-series in interception time period of DAC2F1141056 bearing used for verification
4.2 其他型号轴承的验证
为验证以上规律是否适用于所有汽车轮毂轴承，选取新开发某型汽车轮毂轴承疲劳寿命试验中的振动数据，构建的振动时间序列如图8所示，包含139 318个轴承振动加速度数据，与前述截取规则相同，选取123 148～133 147区间内的10 000个数据，依次分为10个子序列进行分析，计算验证数据10个子序列的形状参数和尺度参数，其变化趋势如图9所示：2种参数的变化规律与表3、图7均具有较好的一致性，即前8个子序列的形状参数和尺度参数在一定范围内波动，第9，10个子序列则发生显著波动，说明本分析方法对汽车轮毂轴承振动时间序列的分析具有通用性。
图8 验证用其他型号汽车轮毂轴承的振动时间序列Fig.8 Vibration time series of automobile hub bearings of other models used for verification
(a) 形状参数 (b) 尺度参数图9 验证用其他型号汽车轮毂轴承截取时间段内不同子序列的形状参数和尺度参数Fig.9 Shape parameters and scale parameters of different sub-series in interception time period of automobile hub bearings of other models used for verification
振动时间序列演变的统计规律对于滚动轴承在线监测装置的应用有良好的工程价值，相较于现有滚动轴承剩余寿命预测方法，本文提出的方法容易进行工程化，且提前预报的时间长度长于现有研究方法。但是，本研究尚在进行中，后续应在以下方面作进一步研究：
1)扩展振动时间序列数据源。从汽车轮毂轴承扩展至其他类型或其他应用领域的轴承，以进一步验证本方法的通用性。
2)判定标准的界定。通过其他理论建立关联性，以确定形状参数和尺度参数变化程度的量化判定标准。
3)实现滚动预测分析。本文为了方便定性分析，采用自下而上的分析方法，即以失效节点倒推截取数据；
实际应用中可采用程序化实现滚动预测，在获取最初1 000个数据后即可计算表征轴承起始状态的形状参数和尺度参数，当测取到1 001个数据时去除第1个数据，采用2～1001个数据计算形状参数和尺度参数表征更新的轴承状态，依次循环，得到形状参数与尺度参数的连续演变规律，最终实现工程化。