时保吉,栗永非
(新乡职业技术学院,河南 新乡 453006)
随着科学技术的快速发展,为确保系统安全可靠运行,现代机械对滚动轴承性能提出了更加严格的要求。其旋转精度、振动、摩擦力矩波动性以及寿命等直接影响主机的可靠性和使用寿命[1-2]。对滚动轴承摩擦力矩变异过程的实时评估与预报,可以为滚动轴承摩擦磨损寿命和可靠性分析提供依据。滚动轴承摩擦力矩的时间序列属于乏信息过程[3-4]。经典统计理论依赖于已知的概率分布与大样本数据,难以有效地解决滚动轴承摩擦力矩不确定性的动态预报问题。
目前,关于滚动轴承的摩擦力矩的研究已取得一些成果。文献[5]提出了最大熵方法,对HTKB轴承摩擦力矩分布特征参数进行了静态评估与预测,并由试验数据验证了模型的正确性;文献[6]以模糊集合理论为基础,从小样本数据入手,建立了航天轴承摩擦力矩参数的经验概率密度函数及其模糊预报模型,并对摩擦力矩的均值上界和最大值上界进行预报,试验研究表明,预报结果可以满足航天工程的要求;文献[7]提出了理论模型和试验方法来表征改进的推力球轴承在混合和全油膜润滑条件下的摩擦力矩,并研究了润滑剂黏度对摩擦力矩的影响;文献[8]提出了转盘轴承启动力矩变异系数的自助最大熵评估方法,但以上成果尚未涉及摩擦力矩不确定性的动态预报。
下文融合了GM(1,1)模型完善的预报机制和自助法的模拟再抽样能力的优点,在乏信息条件下,提出了滚动轴承摩擦力矩不确定性的动态灰自助GBM(1,1)预报方法,通过摩擦力矩不确定性试验,证实了该方法的正确性与有效性。
假设所研究的滚动轴承摩擦力矩为随机变量x。在滚动轴承试验期间,对其摩擦力矩进行定期采样分析,假设获得了该性能的R个时间单元的数据。令Xr为第r个时间单元采集的数据,并构成第r个时间序列Xr
Xr={xr(n)};n=1,2,…,r=1,2,…,R,
(1)
式中:xr(n)为Xr中的第n个原始数据;n为数据序号,即时刻。
借助与时刻n紧邻时刻的前m个数据(包括时刻n的数据)评估时刻n的不确定度。在时刻n只考虑m个数据的原因是:在动态测量中,不确定度随时间而变化,m越大,包含的旧信息越多,估计的不确定度的误差越大。灰色系统理论的灰预测模型GM(1,1)和自助法要求m的最小值为4。
在时刻n,从Xr中抽取m个数据,形成一个小样本数据子序列Xrm
Xrm={xrm(u)+qr};
u=n-m+1,n-m+2,…,n,n≥m,
(2)
式中:u为时间变量;qr为常数,在灰预测GM(1,1)中,若xrm(u)<0,qr应满足xrm(u)+qr≥0;若xrm(u)≥0,则取qr=0。
在时刻n,从Xrm中等概率可放回地随机抽取一个数据,抽取m次,得到一个大小为m的自助样本[9]。连续重复B步,得到B个自助再抽样样本,用向量表示为
YrBootstrap=(Yr1,Yr2,…,Yrb,…,YrB),
(3)
Yrb={yrb(u)};
u=n-m+1,n-m+2,…,n,
n≥m,b=1,2,…,B,
式中:Yrb为第r个时间单元中的第b个自助样本:yrb(u)为Yrb中的第u个自助再抽样数据。
根据灰预测模型GM(1,1)[10],设Yrb的一次累加生成序列向量为
(4)
u=n-m+1,n-m+2,…,n,
n≥m,b=1,2,…,B。
灰生成模型为
(5)
式中:cr1和cr2是待估系数(cr1≠0)。
用增量代替微分,即
xrb(u+1)-xrb(u)=yrb(u+1),
(6)
式中:Δu取单位时间间隔1,设均值生成序列向量为
Zrb={zrb(u)}={0.5xrb(u)+0.5xrb(u-1)};
(7)
u=n-m+2,n-m+3,…,n,
n≥m,b=1,2,…,B。
初始条件xrb(n-m+1)=yrb(n-m+1) 下, (5) 式的最小二乘解为
(8)
D=(-Zrb,I)Τ,I=(1,1,…,1) 。
时刻w=n+1的预测值为
(9)
因此,可以获得第r个时间单元在w=n+1时刻的B个数据,构成如下序列向量
(10)
用统计理论的直方图方法,建立第r个时间单元在时刻w的概率密度函数
frw=frw(xrm) ,
(11)
式中:xrm是描述第r个时间单元中被测数据xrm(w)的随机变量。
假设显著性水平为α,则置信水平为
P=(1-α)×100%,
(12)
且在时刻n,在置信水平P下,xrm的估计区间为
[XrL,XrU]=[XrL(w),XrU(w)]=
[Xrα/2,Xr1-α/2],
(13)
式中:Xrα/2为xrm对应于概率α/2的值;Xr1-α/2为xrm对应于概率1-α/2的值;XrU,XrL分别为估计区间的上、下边界。
定义n时刻的扩展不确定度,即瞬时不确定度为
Ur=Ur(w)=XrU-XrL。
(14)
显然,灰自助预报GBM(1,1)用w=n+1时刻的预报值描述n时刻的瞬时不确定度。定义 (14) 式的不确定度为时刻n的函数,即动态不确定度,其不同于经典统计方法的静态不确定度,而是随着时刻n变化而变化。
在动态测量过程中,假设每个时间单元中,测量数据的总数为nr=N。如果有hr个数据落在估计区间[XrL,XrU]之外,则定义参数PrB为
(15)
式中:PrB是给定置信水平P下,对第r个时间单元估计结果的可靠度,用于描述GBM(1,1)预报的可信度。一般PrB≠P。根据PrB的定义知,PrB越大越好,最好是PrB≥P。
由[XrL,XrU]和PrB的定义可知,在w时刻,P越大,则Ur越大。若P=100%,则Ur可以取最大值。必须注意的是,Ur越大,[XrL,XrU]越偏离真值,进而估计结果越失真,因此,第r个时间单元的平均不确定度Urmean需满足以下条件[11]
(16)
Urmean|m,B,P→min。
(17)
Urmean实际是一个统计量,是动态测量过程中变量不确定度的均值,可以作为动态测量过程中随机变量波动状态的评价指标。第r个时间单元表示随机变量时间历程的第r个阶段。
从 (16)~(17) 式可以看出,最合适的评估结果是在PrB=100%的条件下,Urmean得到最小值。因此,在工程实践中,应结合具体的研究对象,选择满足条件的参数m,B和P。
Urmean对滚动轴承摩擦力矩变异的动态预报有重要作用。Urmean越小,摩擦力矩波动越小,轴承性能变异越轻微;否则,摩擦力矩波动越大,轴承性能变异越严重。据此可以实时预报轴承性能随时间的恶化程度,为及时采取相应措施、消除安全隐患提供决策建议。
试验轴承型号为608,摩擦力矩测量仪型号为M9908A。摩擦力矩主轴的转速为5 r/min,轴承外圈承受的轴向载荷为15 N。取0.005 86 s为一个采样间隔单位,获得了2 048个摩擦力矩数据。取前2 000个数据作为研究对象,将其分为4组,即得到R=4个时间单元,每个时间单元N=500个数据。摩擦力矩的时间序列如图1所示。由图可知,在给定的载荷下,摩擦力矩具有变量不确定性的特点,应对不同时间单元摩擦力矩动态不确定度进行预报,从而实现滚动轴承摩擦力矩不确定性的连续评估,保证系统安全平稳运行。
图1 摩擦力矩的时间序列X1-X4
考虑到模型中m,B和P对Urmean的影响,以第1个时间单元的数据X1为例进行说明。
对乏信息预报而言,B太小,虽然预报速度快且占用内存少,但自助再抽样不充分,计算结果的波动性增大,直接影响预报效果;B太大,降低了预报速度且占用内存过多,不利于在线预报和控制。另外,B过大并不会提高预报的准确性,其有一个极限值。通过大量计算可知,当B=500时,计算结果不再波动且计算时间最短。因此,试验中令B=500,改变m和P,基于灰预测模型GM(1,1)对时间序列X1进行自助预报,结果见表1。
表1 m,B和P的选取对Urmean的影响
由表可知,在P1B=100%,且m=6,P=99%时,U1mean取最小值0.022 54。因此,对试验中其他3个时间单元的试验数据均采取m=6,P=99%,B=500进行分析。
分别计算4个时间单元的平均不确定度,结果如图2所示。对比图1和图2可以看出,对于第r个时间单元,数据波动越剧烈,Urmean的值就越大,摩擦力矩变异越显著。在前2个时间单元中,数据波动很小,Urmean(r=1, 2)在0.02~0.025之间,接近0,摩擦力矩变异最小;而在后2个时间单元中存在局部脉冲,第3个时间单元的脉冲波动最剧烈,U3mean最大,摩擦力矩变异最显著;第4个时间单元的脉冲强度弱于第3个时间单元,U4mean介于U2mean与U3mean之间。显然,第r个时间单元的平均不确定度Urmean较好地实现了摩擦力矩不确定性的实时预报,揭示了摩擦力矩在各阶段时间历程的变异程度。
图2 4个时间序列的平均不确定度
每个时间单元的振动加速度时间序列在各个时刻的估计区间如图3所示。由此获得的时间序列的预报可靠度如图4所示。
图3 时间序列的估计区间
图4 时间序列的预报可靠度
由图4可知,预报可靠度的最大值为100%,最小值为99.79%,表明图3中轴承各个时刻的摩擦力矩数据至少有99.79%包络在估计区间内,证明所建立的预报模型能够有效地预测轴承摩擦力矩的变化范围,为及时控制失效提供了有利条件。
基于GBM(1,1)的滚动轴承摩擦力矩不确定性的动态预报方法,解决了概率分布和变化趋势均未知的条件下摩擦力矩变异的实时预测问题。该方法可实时追踪各种随机变量的波动和趋势项的变化,实现了乏信息条件下动态不确定度的有效预报,补充了现有不确定性评估理论。