, 姜楠, 王立平 清华大学 机械工程系, 精密超精密制造装备及控制北京市重点实验室, 北京 100084
收稿日期：2018-07-20
基金项目：国家科技重大专项（2016ZX04004-004）
作者简介：张云(1974-), 男, 副研究员。E-mail:zhangyun@tsinghua.edu.cn
摘要：针对数控转台性能退化过程较为缓慢且存在波动的特点，采用Wiener过程进行数学建模，利用极大似然法结合试验数据分析模型中各未知参数的特点。根据分析结果将模型中的漂移系数和扩散系数设为恒定值，假设性能退化初值服从正态分布，得出数控转台的相关可靠性函数，并与基于伪寿命分布法得出的结果进行比较。结果表明：基于Wiener过程得出的结果较基于伪寿命分布法得出的结果更加符合实际情况。
关键词：数控转台Wiener过程极大似然估计可靠性估计
Reliability analysis of NC rotary table based on a Wiener process for extremely small samples
ZHANG Yun
, JIANG Nan, WANG LipingBeijing Key Laboratory of Precision/Ultra-Precision Manufacturing Equipments and Control, Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China
Abstract: The performance of numerical control (NC) rotary tables degrades slowly and oscillatorily. A Wiener process model is developed to improve the reliability estimate of NC rotary tables using the maximum likelihood method to analyze the measured characteristics of each parameter with fixed drift and diffusion coefficients assuming that the initial performance degradation values obey a normal distribution. The model gives the correlation reliability function for the NC rotary table. Comparison of the results with those of the pseudo life distribution method shows that the Wiener process model more accurately predicts the actual characteristics.
Key words: numerical control (NC) rotary tableWiener processmaximum likelihood estimationreliability estimate
利用性能退化数据进行产品的可靠性评估是当前较为通用的做法。Nelson[1]、Lu与Meeker等[2-3]在基于性能退化数据的可靠性评估方面做了很多工作；其中使用较多的数学模型有回归模型和Wiener过程模型。对于退化过程存在随机波动的产品，由于回归模型的方差是常量，无法对退化过程中退化量的波动性进行建模，只能对大样本下不同样本的回归方程系数进行建模来评估产品总体的退化失效随机性和可靠性；而Wiener过程数学模型是由Brown运动和线性漂移项合成的，对产品的性能退化过程波动性具有很好的建模效果，因此常用于具有波动性的性能退化过程的建模。
彭宝华等[4]采用Wiener过程对电容器的电容值耗损量建立了数学模型，并对其可靠性进行了评估；朱德馨等[5]将Wiener过程模型中的轴承退化量初值和漂移系数随机化，结合Bootstrap方法对极小样本下高速列车轴承的可靠性进行了建模；王明磊等[6]采用二元Wiener过程对小样本的电主轴可靠性进行了分析；Wang[7]基于Wiener过程对桥横梁的可靠性进行了研究；Tang等[8]采用Wiener过程对发光二极管的亮度建立了模型，并评估其可靠性；Nicolai等[9]利用Wiener过程分析了钢构表面锈迹的扩散；Whitmore[10]利用Wiener过程建立了晶体管增益的退化轨迹模型。
基于Wiener过程的可靠性评估有多种不同的实现方法。如果假设未知参数均为常数，根据所有试验数据构造似然函数来估计未知参数，虽然计算简单但是无法体现不同样本之间的差异；但若假设所有参数均服从正态分布，然后使用所有数据构造似然函数估计未知参数，将会使得评估过程异常复杂。本文先根据样本数据估计退化模型参数，然后根据参数的估计结果假设其恒定或服从正态分布，根据以上假设估计模型参数，最后进行转台可靠性估计。
1 Wiener过程退化模型Wiener退化过程的数学表达式为
| $y\left( {t;v,{\delta ^2}} \right) = a + vt + \delta B\left( t \right).$ | (1) |
1）独立增量性。每个时刻的性能参数的观测值都服从独立同分布。
2）稳定增量和正态性。研究样本属于退化失效型产品，即某个性能参数随着时间的推移表现为单调上升或者下降的趋势，且性能增量服从正态分布。
3）B(0)=0，从而B(ti)~N(0，ti)。
产品的某个性能指标达到某个特定水平时失效，则其失效时间T可以定义为产品的性能退化量X首次达到失效阈值Df所用的时间，即
| $T = \inf \left\{ {t\left| {X\left( t \right) = {D_{\rm{f}}},t \ge 0} \right.} \right\}.$ | (2) |
| $f\left( t \right) = \frac{{{D_{\rm{f}}} - a}}{{\sqrt {2{\rm{ \mathsf{ π} }}{\delta ^2}{t^3}} }} \cdot \exp \left( { - \frac{{{{\left( {{D_{\rm{f}}} - a - vt} \right)}^2}}}{{2{\delta ^2}t}}} \right).$ | (3) |
| $\begin{array}{*{20}{c}}{R\left( t \right) = \mathit{\Phi }\left( {\frac{{{D_{\rm{f}}} - a - vt}}{{\delta \sqrt t }}} \right) - }\\{\exp \left( {\frac{{2v\left( {{D_{\rm{f}}} - a} \right)}}{{{\delta ^2}}}} \right)\mathit{\Phi }\left( {\frac{{ - \left( {{D_{\rm{f}}} - a} \right) - vt}}{{\delta \sqrt t }}} \right).}\end{array}$ | (4) |
| $\Delta {y_i} = v\Delta {t_i} + \delta \Delta {B_i},i = 1,2, \cdots ,N.$ | (5) |
| $\Delta {y_i} \sim N\left( {v\Delta {t_i},{\delta ^2}\Delta {t_i}} \right),i = 1,2, \cdots ,N.$ | (6) |
| $L\left( {v,\delta } \right) = \sum\limits_{i = 1}^N {\frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}{\delta ^2}\Delta {t_i}} }}{{\rm{e}}^{ - \frac{{{{\left( {\Delta {y_i} - v\Delta {t_i}} \right)}^2}}}{{2{\delta ^2}\Delta {t_i}}}}}} .$ | (7) |
| $\begin{array}{*{20}{c}}{\ln L\left( {v,\delta } \right) = }\\{\sum\limits_{i = 1}^N {\left[ { - \ln \sqrt {2{\rm{ \mathsf{ π} }}\Delta {t_i}} - \ln \delta - \frac{{{{\left( {\Delta {y_i} - v\Delta {t_i}} \right)}^2}}}{{2{\delta ^2}\Delta {t_i}}}} \right]} .}\end{array}$ | (8) |
| $\frac{{\partial \ln L}}{{\partial v}} = \sum\limits_{i = 1}^N {\frac{{\Delta {y_i} - v\Delta {t_i}}}{{{\delta ^2}\Delta {t_i}}}} = 0,$ | (9) |
| $\frac{{\partial \ln L}}{{\partial \delta }} = \sum\limits_{i = 1}^N {\left( { - \frac{1}{\delta } + \frac{{{{\left( {\Delta {y_i} - v\Delta {t_i}} \right)}^2}}}{{{\delta ^3}\Delta {t_i}}}} \right)} = 0.$ | (10) |
| $\hat v = \frac{1}{N}\sum\limits_{i = 1}^N {\frac{{\Delta {y_i}}}{{\Delta {t_i}}}} ,$ | (11) |
| $\hat \delta = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {\frac{{{{\left( {\Delta {y_i} - \hat v\Delta {t_i}} \right)}^2}}}{{\Delta {t_i}}}} } .$ | (12) |
1.2 模型参数分布假设假设样本初值服从正态分布，漂移系数和扩散系数为恒定值，即a~N(a0，σ2a)，
当漂移系数和扩散系数恒定、初值a~N(a0，σ2a)时，结合式（4），产品的可靠度函数为
| $\begin{array}{*{20}{c}}{R\left( t \right) = \int {\left( {\mathit{\Phi }\left( {\frac{{{D_{\rm{f}}} - \widehat {{a_0}} - \hat vt}}{{\hat \delta \sqrt t }}} \right) - } \right.} }\\{\exp \left( {\frac{{2\hat v\left( {{D_{\rm{f}}} - \widehat {{a_0}}} \right)}}{{{{\hat \delta }^2}}}} \right),}\\{\mathit{\Phi }\left( {\frac{{ - \left( {{D_{\rm{f}}} - \widehat {{a_0}}} \right) - \hat vt}}{{\hat \delta \sqrt t }}} \right)p\left( a \right){\rm{d}}a.}\end{array}$ | (13) |
概率密度函数为
| $\begin{array}{*{20}{c}}{f\left( t \right) = \int {\left\{ {\frac{{{D_{\rm{f}}} - \widehat {{a_0}}}}{{\sqrt {2{\rm{ \mathsf{ π} }}{{\hat \delta }^2}{t^3}} }} \cdot } \right.} }\\{\left. {\exp \left( { - \frac{{{{\left( {{D_{\rm{f}}} - \widehat {{a_0}} - \hat vt} \right)}^2}}}{{2{{\hat \delta }^2}t}}} \right)} \right)p\left( a \right){\rm{d}}a.}\end{array}$ | (14) |
2 试验验证采用某型号数控转台进行可靠性试验，在不同转台角度θ下得到的转台回转定位精度误差退化量见表 1。回转定位精度误差退化量的退化趋势见图 1。
表 1 数控转台性能退化数据Ⅰ(部分)
| t/h | 精度误差退化量/(″) | |||
| θ＝30° | θ＝45° | θ＝60° | θ＝75° | |
| 50 | 8.5 | 9.0 | 2.0 | 5.0 |
| 100 | 11.5 | 2.0 | 7.0 | 3.5 |
| 150 | 10.5 | 0.5 | 8.5 | 2.5 |
| 200 | 7.0 | 6.0 | 5.5 | 7.0 |
| 250 | 7.5 | 6.0 | 2.5 | 7.5 |
| 300 | 2.0 | 2.0 | 6.5 | 10.0 |
| 350 | 2.5 | 1.5 | 9.0 | 10.0 |
| 400 | 4.0 | 6.0 | 7.5 | 13.5 |
| 450 | 7.0 | 11.0 | 12.0 | 17.0 |
| 500 | 12.0 | 17.0 | 20.5 | 22.5 |
| 550 | 14.5 | 17.0 | 19.0 | 23.0 |
| 600 | 30.5 | 39.0 | 42.0 | 44.5 |
| 650 | 24.0 | 34.0 | 38.5 | 40.5 |
| 700 | 27.0 | 39.5 | 41.5 | 41.5 |
| 800 | 23.0 | 36.5 | 42.0 | 41.5 |
| 850 | 28.0 | 39.0 | 39.0 | 42.5 |
表选项
 |
| 图 1 数控转台回转定位精度误差退化量趋势图 |
| 图选项 |
从图 1可以看出，数控转台的回转精度误差退化量随时间增加不仅有上升趋势也有波动趋势，反映了其性能逐渐退化，且由于试验环境、载荷以及测量误差等因素的影响，在退化过程中性能存在一定的随机波动。这种情况适合采用Wiener过程对其进行建模评估。
由1.2节的计算过程并结合表 1中的数据，得出Wiener模型的各参数，见表 2。
表 2 表 1中各组数据退化量初值、漂移系数和扩散系数
| i | θi/(°) | ai | vi | δi |
| 1 | 30 | 8.5 | 0.026875 | 0.717839 |
| 2 | 45 | 9 | 0.039375 | 0.917122 |
| 3 | 60 | 2 | 0.045938 | 0.893250 |
| 4 | 75 | 5 | 0.046875 | 0.772463 |
表选项
根据表 2数据，计算各个参数的平均值、标准差以及变异系数，见表 3。
表 3 模型参数的平均值、标准差以及变异系数
| a | v | δ | |
| 平均值 | 6.125000 | 0.039766 | 0.825169 |
| 标准差 | 3.275541 | 0.009219 | 0.095548 |
| 变异系数 | 0.534782 | 0.231829 | 0.115792 |
表选项
由表 3可以看出，漂移系数和扩散系数比较稳定，可以将它们设为恒定值，而退化量初值的变异系数较大，可假设它服从正态分布，即a~N(6.125，3.275 5412)，其概率密度函数
| $\begin{array}{*{20}{c}}{R\left( t \right) = 1.002\;71\int_{ - 3.701\;623}^{15.951\;623} {\left( {\mathit{\Phi }\left( {\frac{{{D_{\rm{f}}} - \hat a - \hat vt}}{{\hat \delta \sqrt t }}} \right) - } \right.} }\\{\left. {\exp \left( {\frac{{2\hat v\left( {{D_{\rm{f}}} - \hat a} \right)}}{{{{\hat \delta }^2}}}} \right)\mathit{\Phi }\left( {\frac{{ - \left( {{D_{\rm{f}}} - \hat a} \right) - \hat vt}}{{\hat \delta \sqrt t }}} \right)} \right) \cdot }\\{\frac{1}{{\widehat {{\sigma _a}}}}\varphi \left( {\frac{{a - 6.125}}{{3.275\;541}}} \right){\rm{d}}a.}\end{array}$ | (15) |
 |
| 图 2 数控转台可靠度函数曲线（表 1数据） |
| 图选项 |
根据式（14），概率密度函数为
| $\begin{array}{*{20}{c}}{f\left( t \right) = 1.002\;71\int_{ - 3.701\;623}^{15.951\;623} {\left( {\frac{{{D_{\rm{f}}} - \widehat {{a_0}}}}{{\sqrt {2{\rm{ \mathsf{ π} }}{{\hat \delta }^2}{t^3}} }} \cdot } \right.} }\\{\left. {\exp \left( { - \frac{{{{\left( {{D_{\rm{f}}} - \widehat {{a_0}} - \hat vt} \right)}^2}}}{{2{{\hat \delta }^2}t}}} \right)} \right) \cdot }\\{\frac{1}{{\widehat {{\sigma _a}}}}\varphi \left( {\frac{{a - 6.125}}{{3.275\;541}}} \right){\rm{d}}a.}\end{array}$ | (16) |
 |
| 图 3 数控转台失效概率密度函数曲线（表 1数据） |
| 图选项 |
从图 2和3可以看出，基于Wiener过程的数控转台可靠性估计结果要高于基于伪寿命分布法的可靠性估计结果。在开始的一段时间内，基于Wiener过程的可靠性估计结果显示产品出现失效的概率非常低；而基于伪失效寿命分布法的结果表明产品出现失效的概率较大。从试验过程以及厂家反馈的实际生产使用情况来看，转台在刚开始使用时通常不会出现失效的情况，因此从实际使用结果的角度而言，基于Wiener过程的可靠性评估结果更符合实际情况。从理论上看，数控转台的性能在使用过程中由于环境、试验条件、使用方法以及测量等因素的影响具有随机波动性。对于退化过程存在随机波动的产品，由于回归模型的方差是常量，无法对退化过程中退化量的波动性进行建模，且本文的例子由于属于极小子样情况，也无法对大样本下不同样本的回归方程系数进行建模来评估产品总体的退化失效随机性。因此，基于伪寿命分布法得出的产品可靠性必然小于实际值。Wiener过程数学模型是由Brown运动和线性漂移项合成的，对产品的性能退化过程波动性具有很好的建模效果。因此，从理论分析的角度而言，基于Wiener过程的可靠性评估结果也更符合实际情况。
为了使上一段的论述更具说服力，增加一组试验数据（见表 4），并得到其可靠度函数曲线（见图 4）和失效概率密度函数曲线（见图 5）。图 4和5所示的结果同样表明：基于Wiener过程的数控转台可靠性估计结果要高于基于伪寿命分布法的可靠性估计结果，与上一段的结论一致。可见，对于数控转台的可靠性分析而言，基于Wiener过程的分析结果比基于伪寿命分布法的分析结果更加符合实际情况。
表 4 数控转台性能退化数据Ⅱ(部分)
| t/h | 精度误差退化量/(″) | ||
| θ＝15° | θ＝90° | θ＝105° | |
| 50 | 6.8 | 1.0 | 1.5 |
| 100 | 5.0 | 1.5 | 6.0 |
| 150 | 4.5 | 0.0 | 7.5 |
| 200 | 1.0 | 3.0 | 1.5 |
| 250 | 3.0 | 3.0 | 2.0 |
| 300 | 3.5 | 8.0 | 9.0 |
| 350 | 2.0 | 9.0 | 8.5 |
| 400 | 2.0 | 15.0 | 13.0 |
| 450 | 4.5 | 19.5 | 18.0 |
| 500 | 7.5 | 25.0 | 24.5 |
| 550 | 11.0 | 24.0 | 24.0 |
| 600 | 16.5 | 45.0 | 35.5 |
| 650 | 12.5 | 41.0 | 31.0 |
| 700 | 14.0 | 42.5 | 35.0 |
| 800 | 12.0 | 42.0 | 34.0 |
| 850 | 16.5 | 44.0 | 33.0 |
表选项
 |
| 图 4 数控转台可靠度函数曲线(表 4数据) |
| 图选项 |
 |
| 图 5 数控转台失效概率密度函数曲线(表 4数据) |
| 图选项 |
3 结论本文根据数控转台性能退化的特点，采用Wiener过程对数控转台的性能退化数据进行可靠性建模。利用极大似然法进行模型参数估计，结合试验数据分析结果将退化模型中的性能初值假设为服从正态分布的随机变量而不是定值，充分考虑了数控转台由于各种外界原因产生的个体差异性问题，提高了可靠性评估的精度。最后，将用本文方法得到的可靠性评估结果与基于伪寿命分布法的评估结果进行比较分析，结果表明基于Wiener退化过程的分析结果更加符合实际情况。
参考文献
| [1] | NELSON W. Analysis of performance-degradation data from accelerated tests[J]. IEEE Transactions on Reliability, 1981, R-30(2): 149-155. DOI:10.1109/TR.1981.5221010 |
| [2] | LU C J, MEEKER W Q. Using degradation measures to estimate a time-to-failure distribution[J]. Technometrics, 1993, 35(2): 161-174. DOI:10.1080/00401706.1993.10485038 |
| [3] | MEEKER W Q, ESCOBAR L A, LU C J. Accelerated degradation tests:Modeling and analysis[J]. Technometrics, 1998, 40(2): 89-99. DOI:10.1080/00401706.1998.10485191 |
| [4] | 彭宝华, 周经纶, 金光. 综合多种信息的金属化膜电容器可靠性评估[J]. 强激光与粒子束, 2009, 21(8): 1271-1275. PENG B H, ZHOU J L, JIN G. Reliability assessment of metallized film capacitor using multiple reliability information sources[J]. High Power Laser and Particle Beams, 2009, 21(8): 1271-1275. (in Chinese) |
| [5] | 朱德馨, 刘宏昭. 随机性能退化下极小样本高速列车轴承的可靠性评估[J]. 机械科学与技术, 2013, 32(10): 1499-1504. ZHU D X, LIU H Z. Reliability evaluation of high-speed train bearing based on stochastic performance deterioration with minimum sample[J]. Mechanical Science and Technology for Aerospace Engineering, 2013, 32(10): 1499-1504. (in Chinese) |
| [6] | 王明磊, 原大宁, 刘宏昭. 二元Wiener过程下的小样本电主轴可靠性分析[J]. 机械科学与技术, 2017, 36(2): 279-285. WANG M L, YUAN D N, LIU H Z. Reliability analysis of motorized spindle with small sample based on two-dimensional Wiener process[J]. Mechanical Science and Technology for Aerospace Engineering, 2017, 36(2): 279-285. (in Chinese) |
| [7] | WANG X. Wiener processes with random effects for degradation data[J]. Journal of Multivariate Analysis, 2010, 101(2): 340-351. DOI:10.1016/j.jmva.2008.12.007 |
| [8] | TANG J, SU T S. Estimating failure time distribution and its parameters based on intermediate data from a Wiener degradation model[J]. Naval Research Logistics, 2008, 55(3): 265-276. DOI:10.1002/(ISSN)1520-6750 |
| [9] | NICOLAI R P, DEKKER R, NOORTWIJK J M. A comparison of models for measurable deterioration:An application to coatings on steel structures[J]. Reliability Engineering and System Safety, 2007, 92(12): 1635-1650. DOI:10.1016/j.ress.2006.09.021 |
| [10] | WHITMORE G A. Estimating degradation by a Wiener diffusion process subject to measurement error[J]. Lifetime Data Analysis, 1995, 1(3): 307-319. DOI:10.1007/BF00985762 |
| [11] | 郭琦.基于性能退化数据的可靠性评估方法研究[D].广州: 华南理工大学, 2015. GUO Q. Research on reliability assessment methods based on performance degradation data[D]. Guangzhou: South China University of Technology, 2015. (in Chinese) |
| [12] | 彭宝华.基于Wiener过程的可靠性建模方法研究[D].长沙: 国防科学技术大学, 2010. PENG B H. Research on reliability modeling methods based on Wiener process[D]. Changsha: National University of Defense Technology, 2010. (in Chinese) |
