数控机床是典型的复杂机电液可修系统，可靠性受到工业应用和学术研究的广泛关注^[1-2]。不同于故障即报废的简单机电产品，维修活动会对数控机床可靠性产生影响。在工厂中，数控机床的维修通常以调整或更换一部分零部件为主，故早期研究将维修过程假定为“修复如新”^[3-4]，故障间隔时间服从独立同分布建模。实际上，数控机床的零部件更换调整，可以“修复如新”，但是对于复杂整机而言，维修前后可靠度并没有很大改变。所以，有学者假定维修活动将数控机床的可靠性修复至与故障发生前一刻相同，即“修复如旧”，故障间隔时间服从非齐次Poisson过程^[5-6]。实际数控机床在维修后并不会处于这两种极端，而是处于中间状态即不完全维修状态^[7]。

目前，对于不完全维修状态的机床可靠性研究通常是基于广义更新过程与Kijima模型，文[8-10]等都采用Kijima模型对数控机床进行可靠性建模。但是Kijima模型参数求解复杂，指标求解困难。文[7, 11]提出利用对数线性比例强度模型(log-linear proportional intensity model, LPIM)对数控机床可靠性进行评估。文[7]将LPIM与Weibull分析、非齐次Poisson过程以及广义更新过程等方法对比，说明LPIM的结果优于其他方法。但是，数控机床故障强度通常呈浴盆状^[11-12]，LPIM无法描述存在早期故障和退化故障的故障过程。

本文建立了双重对数比例强度模型(superposed LPIM, S-LPIM) 对数控机床可靠性进行建模。该模型考虑了修复活动对机床可靠性的影响，很好地拟合了数控机床故障率为浴盆状的特性。对一批同品牌同型号服役年龄不同的数控加工中心分为2组采用S-LPIM进行可靠性建模，对模型的时间趋势及修复功效进行检验，验证了方法的有效性。

1 数控机床可靠性建模 1.1 LPIM模型说明

LPIM的故障强度函数为

$ \lambda(t)=\lambda_{0}(t) \mathrm{e}^{\mathit \gamma {m(t)}} . $

(1)

其中：λ₀(t)为基准故障强度函数，γ为修复功效因子，m(t)为平均累计故障数。

在LPIM模型中，每次数控机床的维修活动都会影响其可靠性，即故障强度发生改变，采用γm(t)来反映这种改变。

平均累计故障数可以表示为

$ m(t)=\int_{0}^{t} \lambda(\tau) \mathrm{d} \tau. $

(2)

故可以得到

$ \lambda(t)=\frac{\mathrm{d} m(t)}{\mathrm{d} t}=\lambda_{0}(t) \mathrm{e}^{\mathit \gamma m(t)}. $

(3)

求解式(3)可以得到

$ m(t)=-\frac{1}{\gamma} \ln \left[-\gamma\left(\int \lambda_{0}(t) \mathrm{d} t+C\right)\right]. $

(4)

其中C为常数。给定初始条件m(t)=0，可以得到$C=-\frac{1}{\gamma}$。联立式(3)与(4)可以得到

$ \lambda(t)=\frac{\lambda_{0}(t)}{-\gamma\left(\int \lambda_{0}(t) \mathrm{d} t+C\right)}. $

(5)

1.2 数控机床故障时间S-LPIM建模

文[12-13]指出，数控机床与其他机电类产品一样，故障率呈现典型的浴盆状。早期故障期由于设计的不合理、外购件的缺陷、装配的失误等，具有很高的故障率；但是，通过改进设计或更换零件，数控机床的故障率会急剧下降；最后由于元器件老化、磨损，机床的故障率会再次升高。故基准故障强度函数采用双重对数模型，表示为

$ \lambda_{0}(t)=\mathrm{e}^{a_{1}+b_{1} t}+\mathrm{e}^{a_{2}+b_{2} t} . $

(6)

其中a₁、a₂、b₁、b₂为待定参数。将式(6)代入式(4)和(5)中可以求得

$ \begin{gathered} m(t)=-\frac{1}{\gamma} \ln \left[1+\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)\right], \end{gathered} $

(7)

$ \lambda(t)=\frac{\mathrm{e}^{\left(a_{1}+b_{1} t\right)}+\mathrm{e}^{\left(a_{2}+b_{2} t\right)}}{1+\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)} . $

(8)

1.3 K台数控机床强度函数的参数估计

本文采用最小二乘法对K台同型号服役年龄相同的数控机床的S-LPIM模型参数进行估计。收集的故障时间数据是随机截尾的，假设第k台机床的故障观测时间为[0, T_k]，其中T_k为随机截尾时间。将观测期内采集到的K台机床样本故障数据按照时间顺序排列，得到时刻序列。t_{k, i}表示第k台机床发生第i个故障的时刻，第k台机床故障发生的时刻序列为(t_{k, 1}, t_{k, 2}, …, t_{k, i}, …, t_{k, Nk})，t_{k, Nk}<T_k。总故障次数可以表示为

$ N=\sum\limits_{k=1}^{K} N_{k}. $

(9)

将K台数控机床故障时刻重排，得到总时刻序列为(t₁, t₂, …, t_i, …, t_N)，在t_i时刻，观察到的数控机床故障数计为N(t_i)=i/K。

由于累计故障强度m(t)是非线性函数，非线性最小二乘法问题参数采用Gauss-Newton迭代法求解。对累计故障强度m(t)进行Taylor展开，并忽略高阶项，可以得到

$ \begin{gathered} m(t, \boldsymbol{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}})=m\left(t, \boldsymbol{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}_{0}\right)+\frac{\partial m}{\partial \gamma} \delta \gamma+\frac{\partial m}{\partial a_{1}} \delta a_{1}+ \\ \frac{\partial m}{\partial b_{1}} \delta b_{1}+\frac{\partial m}{\partial a_{2}} \delta a_{2}+\frac{\partial m}{\partial b_{2}} \delta b_{2} . \end{gathered} $

(10)

其中: $\boldsymbol{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}=\left[\begin{array}{lllll}\gamma & a_{1} & b_{1} & a_{2} & b_{2}\end{array}\right]^{\mathrm{T}}$ 为待定参数, $\boldsymbol{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}_{0}$ 为$\boldsymbol{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}$ 初始值。

在$t_{i}$ 处的残差为$e_{i}=m\left(t_{i}, \boldsymbol{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}\right)-m\left(t_{i}, \boldsymbol{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}}_{0}\right)$, 最小二乘法期望残差的平方和最小, 残差的平方和可以表示为

$ \begin{gathered} S=\sum\limits_{i=1}^{N} e_{i}^{2}=\sum\limits_{i=1}^{N}\left\{\frac{\partial m}{\partial \gamma} \delta \gamma+\frac{\partial m}{\partial a_{1}} \delta a_{1}+\frac{\partial m}{\partial b_{1}} \delta b_{1}+\right. \\ \left.\frac{\partial m}{\partial a_{2}} \delta a_{2}+\frac{\partial m}{\partial b_{2}} \delta b_{2}\right\}^{2}=\sum\limits_{i=1}^{N}\left(\boldsymbol{D}^{i} \delta \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\right)^{2}. \end{gathered} $

(11)

其中: $\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}=\left[\begin{array}{lllll}\delta \gamma & \delta a_{1} & \delta b_{1} & \delta a_{2} & \delta b_{2}\end{array}\right]^{\mathrm{T}}$ 表示参数误差,$\boldsymbol{D}^{i}=\left[\frac{\partial m\left(t_{i}\right)}{\partial \gamma} \frac{\partial m\left(t_{i}\right)}{\partial a_{1}} \frac{\partial m\left(t_{i}\right)}{\partial b_{1}} \frac{\partial m\left(t_{i}\right)}{\partial a_{2}} \frac{\partial m\left(t_{i}\right)}{\partial b_{2}}\right] 。$

此时, 残差平方和$S$ 为$\delta \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}$ 的函数, 根据最小二乘法法则, 可求出$\delta \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}$, 进而可以求出$\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}$ 的估计$\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}$。通常, 一次计算结果误差较大, 需要多次迭代, 才可以求出理想的$\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}$。

$D^{i}$ 中的偏导数可由式(7)得到。

$ \begin{gathered} \frac{\partial m}{\partial \gamma}=\frac{1}{\gamma^{2}} \ln \left[1+\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)\right]- \\ \frac{\frac{\mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)}{\gamma\left[1+\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)\right]}, \end{gathered} $

(12)

$ \frac{\partial m}{\partial a_{1}}=-\frac{\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)}{\gamma\left[1+\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)\right]}, $

(13)

$ \begin{aligned} \frac{\partial m}{\partial b_{1}}=\frac{\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}{ }^{2}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}} t \mathrm{e}^{b_{1} t}}{\gamma\left[1+\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)\right]}, \end{aligned} $

(14)

$ \frac{\partial m}{\partial a_{2}}=-\frac{\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)}{\gamma\left[1+\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)\right]}, $

(15)

$ \frac{\partial m}{\partial b_{2}}=\frac{\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}^{2}}\left(1-\mathrm{e}^{b_{2} t}\right)+\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}} t \mathrm{e}^{b_{2} t}}{\gamma\left[1+\frac{\gamma \mathrm{e}^{a_{1}}}{b_{1}}\left(1-\mathrm{e}^{b_{1} t}\right)+\frac{\gamma \mathrm{e}^{a_{2}}}{b_{2}}\left(1-\mathrm{e}^{b_{2} t}\right)\right]} . $

(16)

2 数控机床可靠性指标 2.1 可靠度函数

在t_i时刻，即机床发生i个故障后，数控机床的可靠度函数可以表示为

$ R\left(t \mid t_{i}\right)=\mathrm{e}^{-\int_{t_{i}}^{t+t_{i \lambda}(\tau) d r}} . $

(17)

联立式(2)、(7)与(17)，可以得到

$ \begin{aligned} R\left(t \mid t_{i}\right)=& \exp \left\{-\frac{\mathrm{e}^{a_{1}+i \gamma}}{b_{1}}\left[\mathrm{e}^{b_{1}\left(t_{i}+t\right)}-\mathrm{e}^{b_{1} t_{i}}\right]-\right.\\ &\left.\frac{\mathrm{e}^{a_{2}+i \gamma}}{b_{2}}\left[\mathrm{e}^{b_{2}\left(t_{i}+t\right)}-\mathrm{e}^{b_{2} t_{i}}\right]\right\} . \end{aligned} $

(18)

2.2 平均故障间隔时间

故障强度函数λ(t)表示单位时间内故障的次数，则其倒数即可表示一次故障所经过的时间，即瞬时平均故障间隔时间(instantaneous mean time between failure, IMTBF)可表示为

$ \mathrm{IMTBF}=\frac{1}{\lambda(t)} . $

(19)

累计平均故障间隔时间(cumulative mean time between failures, CMTBF)可以表示为

$ \mathrm{CMTBF}=\frac{t_{2}-t_{1}}{\int_{t_{1}}^{t_{2}} \lambda(\tau) \mathrm{d} \tau} . $

(20)

3 故障趋势及修复功效检验

对得到的数控机床可靠性模型的有效性进行验证，采用似然比检验法进行早期故障趋势、退化故障趋势及修复功效检验。联立式(7)与(18)，可以得到第k台数控机床发生第i次故障的条件概率为

$ \begin{aligned} f\left(t_{i} \mid t_{i-1}\right)=\left(\mathrm{e}^{a_{1}+b_{1} t_{i}+\gamma(i-1)}+\mathrm{e}^{a_{2}+b_{2} t_{i}+\gamma(i-1)}\right).\\ \mathrm{e}^{-\frac{\mathrm{e}^{ \mathrm{a}_{1}+(i-1) \gamma}}{b_{1}}\left[\mathrm{e}^{b_{1} t_{i}}-\mathrm{e}^{b{1} t_{i-1}}\right]-\frac{\mathrm{e}^{a_{2}+(i-1) \gamma}}{b_{2}}\left[\mathrm{e}^{b_{2} t_{i}}-\mathrm{e}^{b_{2} t_{i-1}}\right]} . \end{aligned} $

(21)

得到数控机床故障数据的似然函数为

$ L=\prod\limits_{i}^{N} f\left(t_{i} \mid t_{i-1}\right) . $

(22)

3.1 早期故障趋势检验

原假设H₀：机床故障没有显著的早期故障趋势(b₁=0)。备选假设H₁：数控机床在服役早期故障强度存在显著的增强或减弱趋势(b₁≠0)。统计检验量为

$ \varLambda_{1}=-2 \ln \left(\frac{L\left(b_{1}=0 ; \gamma, b_{2}, a_{1}, a_{2}\right)}{L\left(\gamma, a_{1}, b_{1}, a_{2}, b_{2}\right)}\right) \sim \chi_{1}^{2}. $

(23)

3.2 退化故障趋势检验

原假设H₀：机床故障没有显著的退化故障趋势(b₂=0)。备选假设H₁：数控机床在服役后期故障强度存在显著的增强或减弱趋势(b₂≠0)。统计检验量为

$ \varLambda_{2}=-2 \ln \left(\frac{L\left(b_{2}=0 ; \gamma, a_{1}, b_{1}, a_{2}\right)}{L\left(\gamma, a_{1}, b_{1}, a_{2}, b_{2}\right)}\right) \sim \chi_{1}^{2}. $

(24)

3.3 修复功效检验

原假设H₀：维修活动对机床可靠性无显著影响(γ=0)。备选假设H₁：维修活动会改善或恶化机床的可靠性(γ≠0)。统计检验量为

$ \varLambda_{3}=-2 \ln \left(\frac{L\left(\gamma=0 ; a_{1}, b_{1}, a_{2}, b_{2}\right)}{L\left(\gamma, a_{1}, b_{1}, a_{2}, b_{2}\right)}\right) \sim \chi_{1}^{2}. $

(25)

4 实例分析 4.1 模型求解

本文对一段时间内某品牌同型号3台数控加工中心故障数据进行了收集，截尾时间随机。总观察时间为T=13 728 h，收集到故障时间如表 1所示。从使用情况看，第1台数控加工中心服役较早，故障率较高，第2、3台数控加工中心同时服役，故障状况接近，故将它们分为A组和B组进行分析。A组为第1台加工中心，B组包括第2、3台加工中心。对B组2台加工中心进行故障时刻序列重排，得到相应的总故障时刻序列，将2组数控加工中心的故障数据带入模型，并采用最小二乘法估计模型参数，可以得到参数估计值如表 2所示。

2组加工中心因为服役时间不同，所以基准故障强度函数参数a₁、b₁、a₂、b₂差异较大。但因为该用户企业对加工中心的维修保养策略相同，故得到的维修功效因子接近。2组机床的修复因子均为γ < 0，即该用户企业对数控加工中心的维修活动可以提升机床的可靠性水平。

累计故障数的拟合曲线如图 1、2所示，累计故障数实际值与拟合值的相关系数可以表示为

$ r=\sqrt{1-\frac{\sum\limits_{i=1}^{\hat{N}}\left(N\left(t_{i}\right)-m\left(t_{i}\right)\right)^{2}}{\sum\limits_{i=1}^{\hat{N}}\left(N\left(t_{i}\right)-\bar{N}\left(t_{i}\right)\right)^{2}}}. $

(26)

其中：N(t_i)表示t_i时刻观察到的指定加工中心的故障次数，N(t_i)表示累计故障次数的平均值。计算可得，A、B两组加工中心累计故障数实际值与拟合值的相关系数分别为r_A=0.998 0、r_B=0.987 3。2组累计故障数实际值与拟合值的相关系数都接近1，拟合程度较好，表明本文建立的基于S-LPIM的可靠性模型可以恰当描述加工中心的故障过程，相应模型参数的估计值可信度较高。

4.2 可靠性指标分析

由表 2中模型参数的估计值，可以得到两组加工中心的故障强度函数，如图 3所示，2组加工中心的故障强度函数都呈现简单的浴盆状。被观察的加工中心是已经服役一段时间，基本处于快要进入耗损故障期阶段，故障率在后期出现明显升高。

最后一次故障发生后数控加工中心的可靠度是用户企业比较关心的，故对2组数控机加工中心在最后一次故障发生后的可靠度进行了求解，如图 4所示。因为A组加工中心在观察后期故障频发，最后一次故障发生400 h后，可靠度基本都降为0；B组加工中心在最后一次故障发生1 800 h后，可靠度基本降为0。这也验证了A组数控机床因为服役时间较长，故障较高的实验结果。

平均故障间隔时间是数控加工中心重要的可靠性指标，根据2.2节中定义，本文求解了2组数控加工中心的IMTBF，如图 5所示；从观察开始至当前时间的CMTBF，如图 6所示。A、B两组加工中心因在投入使用前，进行了大型的维修保养，所以IMTBF明显地先升高后降低，CMTBF也存在类似的趋势。B组加工中心的IMTBF与CMTBF，在观察后期出现明显下降。两组相比而言，A组的CMTBF低于B组，在相同观察时间，A组机床的故障更多。

4.3 故障趋势与修复功效检验

1) 早期故障趋势检验。2组加工中心早期故障趋势检验值大于显著水平为5%的临界值，即具有足够的理由拒绝原假设，表明机床在服役早期故障强度存在显著的故障增强或降低趋势。

2) 退化故障趋势检验。2组加工中心统计检验量远大于临界值，加工中心可靠度在服役后期故障强度具有显著增强或退化趋势。由2组加工中心S-LPIM模型中b>0及IMTBF与CMTBF分析可以得出，2组加工中心在服役后期可靠性急剧退化。

3) 修复功效检验。2组加工中心统计检验量远大于临界值，表明维修活动对数控加工中心可靠性具有显著影响，维修功效因子小于0，表明维修活动对加工中心的可靠性提升具有积极作用，即会改善数控加工中心的可靠性。

5 结论

本文提出了双重对数线性比例强度模型用于数控机床可靠性模型，可以反映数控机床故障率呈现浴盆状的特性，并定量反映出修复活动对数控机床可靠性的影响。通过对某企业3台数控加工中心分为2组进行分析表明，所分析的A组、B组加工中心的平均故障间隔工作时间都较低，且在观察后期，2组加工中心的可靠度、瞬时故障间隔工作时间都急剧下降，说明被观察的加工中心已经接近耗损故障期。查看加工中心的故障原因，也多为损耗型故障，与从模型得到的结论一致。不同的维修活动对机床可靠性的影响是不同的，本文中并没有区分何种维修活动，仅对维修次数进行了统计，认为维修活动存在统一的修复功效因子。维修活动的类型对数控机床可靠性的影响及相关建模将有待进一步研究。