考虑外界影响的元动作单元故障概率模型

鞠萍华，谷豪东，许文林，黄广全

（重庆大学 机械工程学院，重庆 400044）

对于数控机床而言，其功能是通过各个部件运动来实现的，而部件的运动则是由一个个基本动作完成的.元动作单元结构化分解就是根据这一特点，将产品按照“功能（function）－运动（motion）－动作（action），FMA”进行分解，其中心思想就是将实现产品功能的复杂的综合运动拆分为一系列基本的动作.将分解得到的最基本的动作称为元动作，将实现元动作的单元称为元动作单元（Meta-action unit，MAU）.一直以来，关于数控机床可靠性的研究比较多，但元动作单元可靠性的研究则比较缺失.元动作单元作为组成数控机床的基本构件，研究其故障概率模型将对数控机床组成构件的可靠性研究有着重要作用.

为得到元动作单元故障概率模型，首先需要对元动作单元产生的故障进行分类.在本文研究之前，对元动作单元故障进行分类的方法有很多，这些方法分别根据不同的标准进行分类.有根据元动作单元故障发生的位置，将元动作单元的故障分为内部故障、连接型故障和控制系统故障[1]；还有按照故障发生和发展的进程将故障分为突发性故障和渐发性故障；除此之外，还有根据观察到的故障现象将故障分为动态异常和静态异常两种类型[2].这些分类方式可以由故障结果直接进行判定，分类比较方便，但不利于得到元动作单元的故障概率模型.

除此之外，对各种机电设备进行分类并得到故障概率模型的研究也有很多，这些概率模型分别使用了不同的分布.其中最简单和常用的是选用泊松分布对故障数据进行描述[3-4]，使用泊松分布描述设备故障是最简单直接的.但泊松分布展现的故障率是恒定的，而机电设备则会由于老化导致故障率不断增加，因此单纯使用泊松分布不能对故障数据进行描述.Xie 等[5]提出用指数分布函数对故障率进行描述.然而，这种方法却将老化故障当作随机故障，无法解释老化故障的故障率随时间增加的趋势.Li[6]提出了使用正态分布和威布尔分布.这种方法可以得到故障率分布，但是这种方法却有一个缺点，其认为随机故障的故障率随时间改变.

上述故障概率模型在描述元动作单元故障数据时存在各种各样的缺陷，本文针对元动作单元的特点，建立元动作单元的故障概率模型.

1 元动作单元

元动作是数控机床进行FMA 结构化分解中得到的最基本的动作，其具体的分解步骤如图1 所示.根据元动作的特点，在此将元动作的概念定义为：元动作为机电产品中传递运动和动力的最基本形式，它是机械产品中最小的运动[7].

图1 FMA 结构化分解示意图Fig.1 FMA structured decomposition diagram

单个零件是无法完成一个规定的元动作，需要一个完整的单元才能实现.元动作单元为能够保证机械产品的元动作得以正常运行的所有零件按照装配关系组成的统一整体，且元动作单元在结构上不可再分也无需再分[7].

复杂的机电产品的功能是由部件的运动来实现，部件的运动是由一系列元动作来实现，而完成一个元动作则需要一个元动作单元.因此当产品发生故障时，最终体现在元动作单元上.而元动作单元的故障是指部分动作单元功能失效，造成不能顺利完成规定动作要求的状态[1].当元动作单元出现故障时，会导致部件的运动出现障碍，最终体现在机电产品出现故障.以数控机床为例，数控机床的运动由机床的进给、工作台的回转、刀架的回转、从刀库抓取刀具、主轴的回转等功能实现.而刀架的回转分为刀架的移动和转动，而对于这两个运动来说，需要两个元动作单元来分别实现.当其中一个元动作单元发生故障时，刀具回转将无法实现，进而数控机床无法工作[8].因此研究单个元动作单元的故障概率模型对于研究数控机床的可靠性来说，十分重要.其故障概率模型的建立过程如图2 所示.

图2 元动作单元故障概率模型流程图Fig.2 Flow chart of failure probability model of meta-action unit

2 元动作单元的故障

对于数控机床而言，从安装完成、投入运营到报废的整个过程的故障率曲线应该如图3 所示，呈现出浴盆曲线的形状[9].早期故障的发生原因在于设计和制造阶段的缺陷，并且故障率随着时间慢慢下降；一段时间后，进入偶然故障期，数控机床工作的大部分时间处于偶然故障期，故障率逐渐稳定，在此期间，由于疲劳和磨损的原因，故障率随时间略微增加，并不为常数；随着时间的增加，由于老化等原因，故障率不断增加，数控机床进入耗损故障期[10].数控机床在交付时会经历充分的试验安装，早期故障期较短，其故障率主要集中在浴盆曲线的偶然故障期和耗损故障期.元动作单元是组成数控机床的最基本单元，其故障率曲线也应该集中在偶然故障期和耗损故障期.

图3 浴盆曲线示意图Fig.3 Schematic diagram of bathtub curve

为了确定元动作单元的故障概率模型，需要对其发生的故障进行分类.本文根据元动作单元故障发生的原因将故障分为老化故障和随机故障.元动作单元在运行过程中由于老化失效而造成的故障称为老化故障；由于一些偶然因素而造成的故障则称为随机故障.

元动作单元是一个关于结构的概念，包含的结构种类非常多，发生的不同类型故障也很多，除此之外，引发同一种故障背后的原因也可能有多种，因此根据故障发生的原因将故障进行分类在操作上会有一些复杂.例如，田致富在论文中对刀架转动异响故障树分析，其结构如图4 所示，引起异响的原因有紧固件松动、啮合区有异物和润滑不充分[11].其中啮合区有异物和润滑不充分为偶然因素造成的故障，紧固件松动有可能是因为老化因素造成的故障，也有可能由于操作过程中偶然因素造成的故障.

图4 刀架转动异响故障树Fig.4 Fault tree for abnormal noise of machine tool holder rotation

在对元动作单元故障实际进行分类中，可参考以下方法，老化故障是元动作单元在工作过程中组成零件的技术指标逐渐下降，最终超出允许的范围而发生故障，可以在事先发现征兆，具有可预测性；随机故障是元动作单元由于偶然因素突发的故障，一般没有任何征兆，具有不可预测性.

为对元动作单元发生的故障进行准确分类，从而确定元动作单元故障概率模型，故障数据的收集不仅要记录故障发生的时间，还要记录发生的部位和发生的具体情形，以便于确定引起故障原因，从而对元动作单元发生的故障进行分类.

3 建立元动作单元的故障概率模型

由于两种故障不同的数据分布特点，所以应选择两个不同的概率分布函数来分别对其进行描述.除此之外，两种故障类型分属于两个独立的部分，所以最后得到的用来拟合元动作单元故障概率模型应为两个不同故障概率分布直接相加.

当故障数据具有三个特点时，可以选用泊松分布来进行描述[12]：

①物体在某段时间内发生n 次故障的概率与时间起点无关，仅与该段时间的长短有关.

②在两段相互不重叠的时间内，物体发生故障的次数n1和n2是相互独立的.

③在较短的时间内，物体发生两次或更多故障的概率很小.

元动作单元的随机故障是由元动作单元内外的偶然因素导致的，例如，带传动的带突然断裂，异物进入元动作单元内等.这些因素随机发生，导致元动作单元发生故障，相互之间没有关联.因此发生n 次故障的概率与时间起点无关，仅与时间长度有关；且互不重叠的时间内，发生故障的次数相互独立；另外，元动作单元是针对机械产品的概念，除特殊情况外，机械产品的寿命较长，可靠性也较高，发生故障的概率较低.综上所述，元动作单元随机故障数据的特点符合泊松分布的要求，故随机故障选用泊松分布进行描述.泊松分布的故障率的计算如式（1）所示：

式中：n 为故障发生的次数，λ0表示故障率，t 为时间.

老化故障是由于元动作单元运行过程中老化失效而发生的故障，例如：齿轮磨损，紧固件老化松动等.老化故障的故障率随时间的增长而逐渐增加，因此应采用当参数变化时，概率分布函数形状随之改变的威布尔分布函数.

威布尔分布的分布函数和故障率分别如式（2）和（3）所示：

其中：η 是比例参数，m 是形状参数，t 为时间

当威布尔分布形状参数改变时故障率随时间变化的情况如图5 所示.当m=1，则威布尔分布为指数分布，故障率为常数；当m<1，则故障概率随时间减小；当m>1，则故障率随时间增加[13].故当m>1时，威布尔分布适合描述老化过程.

图5 形状参数不同时威布尔分布故障概率图Fig.5 Failure probability diagram of Weibull distribution with different shape parameters

由于两种类型的故障相互独立，所以元动作单元故障率λ（t）应为两个不同分布的故障率之和，如式（5）所示：

当t=0 时，元动作单元故障率为λ0；当选择m>1 时，随着时间t 逐渐增加，元动作单元故障率随着时间慢慢增加，符合实际情况.

4 元动作单元故障模型参数估计

元动作单元故障概率模型由独立的两个分布组成，泊松分布和威布尔分布，两种分布属于不同的类型，进行参数估计时应将两个分布分别进行估计.由泊松分布的性质可知，其故障率可以由式（6）计算：

式中：n 为元动作单元发生随机故障的次数，ti为发生随机故障的元动作单元的运行时间，tj为未发生随机故障的元动作单元的运行时间，Z1、Z2分别代表两种类型的元动作单元.

常用的威布尔分布参数估计方法有三种，极大似然估计、最小二乘法、矩估计法[14].三种参数估计方法中，矩估计法和极大似然估计法精度高于最小二乘法[15]，本文选用矩估计法进行参数估计.

由威布尔分布的性质可知，其均值M 和方差S2分别由比例参数m 和形状参数η 来决定：

比例参数m 通过式（7）和式（8）消除后得到只含有形状参数η 的方程（9）.

由统计得到的故障数据分别计算一阶样本矩M和二阶样本矩S2得：

对于Γ 函数而言，存在以下性质：

式中：B 函数为Beta 函数.

因此，对式（9）中包含的Γ 函数进行转化可以得到：

对于一组统计得到的故障数据来说，S2和M 确定，方程（9）的右边相当于常数.方程（10）为关于η的一元方程，方程左边为关于η 的函数，则方程的解应为函数f（η）与常数交点的横坐标.由B 函数的性质可知，函数f（η）在（0，+∞）内单调递减，随着η 增加，函数无限接近于0，因此方程（9）有且只有一个解[16].方程可以通过二分法进行数值计算得到唯一确定的形状参数η，求得η 后，根据式（7）则可以求得比例参数m.

5 元动作单元工作负载参数和工作环境参数的估计

元动作单元在工作过程中，其故障率势必会受到外界因素的影响.本文中，故障类型被分为两类，老化故障和随机故障.在外界影响因素中，工作负载和工作环境会对两种故障类型的故障率分别造成影响.为比较各个工作地区的工作负载和工作环境对元动作单元故障率的影响，本文分别采用工作负载参数和工作环境参数.

元动作单元工作过程中，工作负载越大，则老化速度就越快，所以应在威布尔分布中添加工作负载参数.而工作环境会对随机故障的故障率造成影响，应对泊松分布添加工作环境参数[17].

当平均负载不断增大时，设备老化的速度不断加快，比例参数η 往往与工作负载的大小有关，负载越大，比例参数越小[12]，因此工作负载参数Rl应表现在威布尔分布的比例参数上；由于工作环境会直接影响随机故障的故障率，且泊松分布只有故障率λ一个参数，所以工作环境参数Re应直接表现在泊松分布的故障率上.

工作负载参数Rl和工作环境参数Re的估计方法和步骤如图6 所示.当工作负载参数Rl大于1，则说明对该元动作单元来说工作负载小于平均水平，如果小于1 则说明工作负载大于平均水平；且工作负载参数越大，则说明工作负载越小；反之，参数越小，则说明工作负载越大.

图6 确定工作负载参数和工作环境参数流程图Fig.6 Flow chart for determining working load parameters and working environment parameter

工作环境参数Re则相反，如果大于1，则说明工作环境比平均水平更恶劣，如果小于1，则说明工作环境比平均水平更优；除此之外，工作环境参数越大，则说明工作环境更差，反之，则说明工作环境更好.

6 算例分析

某批次数控机床CJK1630 共计116 台平均分成四组，交付用户使用便开始进行现场试验，选取其中重要的转动元动作单元（主轴轴系结构），待所有机床的该元动作单元都出现故障，则数据收集结束，主轴的故障类型如图7 所示.

图7 主轴故障图Fig.7 Spindle failure diagram

将统计得到的故障数据按照故障发生的原因进行分类后，得到的随机故障数据和老化故障数据分别如表1、表2 所示.

表1 元动作单元随机故障数据Tab.1 Failure data for random failures of meta-action units

表2 元动作单元老化故障数据Tab.2 Failure data for aging failures of meta-action units

在本文中，将元动作单元用字母k 进行表示，由于在前文中将元动作单元平均分为四组，将四组元动作单元分别用1，2，3，4 进行编号，且在接下来的表格中均用此进行表示.

如前文所述，随机故障数据用泊松分布进行拟合，老化故障数据用威布尔分布进行拟合.假设同一批次生产的数控机床性能、工作环境和工作负载均相同，那么应该用相同的参数集A=（λ0、η0、m0）对所产生的全部故障数据进行拟合.

对泊松分布进行参数估计，在本例中，没有未发生故障的元动作单元，故tj为0：

采用矩估计法对威布尔分布进行参数估计时，为保证估计的结果精确，对式（9）采用二分法进行参数估计时，方程两端的差值应小于0.001.

对总的老化故障数据最终估计的结果为η0=2.35，m0=2 229.674

绘制老化故障数据的概率密度函数曲线f（t）和频率直方图，结果如图8 所示.

图8 老化故障数据概率密度函数曲线和频率直方图Fig.8 Probability density function curve and frequency histogram of aging failure data

由图8 可知，老化故障数据的概率密度函数曲线和频率直方图具有较好的拟合效果，说明威布尔分布适合用来描述老化故障数据，且提出的故障概率模型有效.

由于生产过程中各种偶然因素的影响，且所处的工作环境和工作负载各不相同，因此，对四组用户的数控机床故障数据应该采用不同参数的概率模型来进行拟合[18].

对泊松分布进行参数估计，tj亦为0：

最终计算的结果如表3 所示.

表3 各组元动作单元老化故障率Tab.3 Aging failure rate of each meta-action unit

采用矩估计法对各组元动作单元故障数据进行威布尔分布参数估计，其结果如表4 所示.

表4 各组元动作单元威布尔分布参数表Tab.4 Weibull distribution parameter table of each meta-action unit

分别绘制4 组元动作单元的老化故障数据概率密度函数曲线f（t）和频率直方图，结果如图9 所示.

由图9 可知，各组元动作单元的老化故障数据概率密度函数曲线和频率直方图都具有较好的拟合效果，说明各组元动作单元的老化故障数据同样适合用威布尔分布进行描述，进一步说明了提出的故障概率模型有效.

图9 各组元动作单元老化故障数据概率密度函数曲线和频率直方图Fig.9 Probability density function curve and frequency histogram of aging failure data of each meta-action unit

工作负载和工作环境影响元动作单元的故障率，对于这两个参数的估计，按图5 所述，经计算，四组外界影响参数如表5 所示.

由表5 中的外界影响参数可知，第1 组和第2组元动作单元，工作负载小于平均水平，工作环境比平均水平更恶劣；第3 组和第4 组元动作单元则相反，工作负载大于平均水平，工作环境比平均水平更优.

表5 外界影响参数Rl 和Re 表Tab.5 External influence parameters Rl and Re

7 结论

1）本文提出了一种由威布尔分布和泊松分布组成的元动作单元故障概率模型，并考虑了工作负载和工作环境对元动作单元故障率的影响，最终发现元动作单元的故障率随时间变化的形式类似于失去早期阶段的浴盆曲线.

2）元动作单元是数控机床完成一个动作的最基本单元，其故障概率模型不仅可以对复杂机电设备的重要运动构件的可靠性分析提供帮助，还可以为运动构件维修保养计划的制定提供依据.

3）本文提出的元动作单元故障率模型对故障数据的处理简单，经算例分析验证有效，并且上述的模型不仅局限于元动作单元，在可靠性领域中，故障来源可以分为老化故障和随机故障的结构或系统均可以使用.