一种基于隐马尔科夫链的数控装备可靠性分析方法

摘要

本发明提供了一种基于隐马尔科夫链的数控装备可靠性分析方法，具体为：1.监测数控装备的动态性能信号，抽取反映数控装备可靠性变动的性能特征参数值；2.构建性能特征参数值的预测模型；3.使用预测模型预测矢量化所需时间点内的性能特征参数值，采用信源编码方法对这些性能特征参数预测值进行矢量化；4.采用离散隐马尔科夫链模型识别数控装备的状态变迁概率；5.利用切普曼－柯尔莫哥洛夫微分方程建立运行状态与状态变迁概率的关系式，推断数控装备处于不同运行状态下的概率，即获得数控装备的可靠性变动规律。本发明能够在数控装备发生故障前准确地分析、评估和预测数控装备可靠性变动情况，避免数控装备故障的发生，提高数控装备的运行可靠性。

说明书

技术领域

本发明所属于数控装备技术领域，具体是一种基于离散隐马尔科夫链模型的数控装备可靠性变动分析方法。

背景技术

可靠性是数控装备的一种内在属性，是指数控装备在规定任务及其条件下一个规定时间内执行其要求任务的各种关键功能的能力。它反映了数控装备使用过程中的质量特性，考虑了数控装备运行时的时间因素和环境因素。

目前，数控装备可靠性分析的一般思路是通过采集和分析数控装备组成零部件和整机系统的设计与制造信息、现场故障数据以及同类产品的历史故障与寿命数据，从数控装备的故障诊断分析入手，根据数控装备工作和功能特点，寻找故障模式和原因，评估数控装备可靠性水平，并提出可靠性改进措施。其中，现场故障数据采集的途径多是长时间跟踪几十台数控装备，故障数据来源于数控装备的现场记录和维修报告。由此可见，现有的数控装备可靠性分析方法需要有长时间工作、数量足够多的同型号数控装备，且要有大量故障与寿命数据；通过分析故障数据，并基于一定的统计判断准则，识别装备的失效概率分布，进而采用数理统计方法推断出数控装备可靠性，找出数控装备使用过程中的薄弱环节。

然而，这种基于故障数据的数控装备可靠性研究面临的主要问题有：一方面，它需要大量的故障与寿命数据，以发现数控装备的寿命分布曲线。实际工程中几乎不可能得到大量的数控装备故障与寿命数据，特别是整机系统的数据。此外，数控装备经过长期运行出现故障，此时装备可能已处于寿命末期，导致所得的可靠性结果缺乏实际指导意义。另一方面，它仅仅考虑了数控装备的正常和故障两种状态，没有考虑数控装备使用过程中性能状态变化特征。例如，数控装备发生故障前往往经历了一系列的劣化状态，具体表现为输出性能参数逐步偏离其规定值，直至最终发生功能失效。这导致用户无法了解数控装备可靠性的实时变动情况。

因此，迫切需要设计与开发一种用于数控装备可靠性领域的新方法，它不需要等待数控装备经过长时间运行而出现故障，就能准确地分析、评估与预测数控装备可靠性变动情况，从而指导用户提前采取合理的预防性措施，防止故障发生。

发明内容

本发明的目的是在于针对现有技术的不足，提供一种基于离散隐马尔科夫链模型的数控装备可靠性变动分析方法，在数控装备发生故障之前准确地分析、评估和预测数控装备可靠性变动情况，找出数控装备使用过程中的薄弱环节等，避免数控装备故障的发生，提高数控装备的运行可靠性。

一种基于隐马尔科夫链的数控装备可靠性分析方法，包括以下步骤：

(1)监测数控装备运行时的动态性能信号，在一个以上的时间点抽取一个以上的反映数控装备可靠性变动的性能特征参数值；

(2)依据抽取的性能特征参数值，构建性能特征参数值的预测模型；

(3)使用预测模型预测T₀个时间点的各个性能特征参数值，对这些预测值进行矢量量化，得到性能特征参数的输出观测序列；

(4)依据输出观测序列，采用离散隐马尔科夫链模型确定数控装备的运行状态变迁概率矩阵；

(5)利用切普曼-柯尔莫哥洛夫微分方程建立运行状态与状态变迁概率矩阵的关系式，推断数控装备处于不同运行状态下的概率。

所述步骤(2)采用非平稳自回归积分滑动平均算法建立第j个性能特征参数值的预测模型$y_{\mathrm{tj}}=\underset{z=1}{\overset{p}{\Sigma }}{\phi }_z{y}_{(t-z)j}+{ϵ}_{\mathrm{tj}}-\underset{z=1}{\overset{q}{\Sigma }}{\phi }_z{ϵ}_{(t-z)j},$其中，t表示预测时间点，p为自回归阶次，q为滑动平均阶次，φ_z为滑动平均系数，测量误差${ϵ}_{\mathrm{tj}}~N(0,{\sigma }_{ϵ}^2).$

所述步骤(3)具体为：性能特征参数的输出观测序列表示为$O^j=\{o_1^j,···,o_{T_0}^j\},$$o_{\tilde{t}}^j\in (v_1^j,,v_N^j),$$\tilde{t}=1,···,T_0,$j＝1，…，m，m为性能特征参数个数，t表示预测时间点，其中，

$v_1^j=1,v_N^j=N,$N为离散区域数；

$v_{i_0}^j=\mathrm{Index}\left(y_{t_{1}j}\right)=\left(\begin{array}{c}1,y_{t_{1}j}\le \mathrm{partition}(1,j)\\ i_0,\mathrm{partition}(i_0-1,j)\le y_{t_{1}j}\le \mathrm{partition}(i_0,j);\\ N,\mathrm{partition}(N-1,j)\le y_{t_{1}j}\end{array}\right)$

i₀＝2，…，N-1，t₁∈(1，…，T₀)，表示第j个性能特征参数在t₁时间点的预测值，$\mathrm{partition}(i_1,j)=h_{i_{1}j},$i₁＝1，…，N-1，表示第j个性能特征参数处于第i₁离散区域的界定值。

所述步骤(4)具体为：令数控装备包含0～k+1种运行状态，依据步骤(3)得到的输出观测序列以及预定义的初始状态变迁概率矩阵A₀、初始观测值概率矩阵B₀和初始概率分布向量π₀，利用鲍姆-韦尔奇(Baum-Welch)算法迭代计算得到数控装备在t时刻的状态变迁模型λ＝(A，B，π)，π表示概率分布向量，

状态变迁概率矩阵$A={\left[a_{\mathrm{cd}}\right]}_{(k+1)\times (k+1)}=\left(\begin{array}{ccccc}{a}_{00}& a_{01}& ···& a_{0k}& a_{0(k+1)}\\ 0& a_{11}& ···& a_{1k}& a_{1(k+1)}\\ ···& ···& ···& ···& ···\\ 0& 0& ···& a_{\mathrm{kk}}& a_{k(k+1)}\\ 0& 0& ···& 0& a_{(k+1)(k+1)}\end{array}\right),$a_cd为数控装备在t时刻处于状态{c}而在t+1时刻处于状态{d}的概率，0≤a_cd≤1，且$\underset{d=0}{\overset{k+1}{\Sigma }}{a}_{\mathrm{cd}}=1;$

观测值概率矩阵$B={\left[b_{\mathrm{cj}}\right]}_{(k+1)\times m}=\left(\begin{array}{cccc}{b}_{01}& b_{02}& ···& b_{0m}\\ b_{11}& b_{12}& ···& b_{1m}\\ ···& ···& ···& ···\\ b_{k1}& b_{k2}& ···& b_{\mathrm{km}}\\ b_{(k+1)1}& b_{(k+1)2}& ···& b_{(k+1)m}\end{array}\right),$b_cj为数控装备在t时刻处于状态{c}时出现观测值$\left\{v_c^j\right\},c=0,···,k+1$的概率，0≤b_cj≤1，且$\underset{j=1}{\overset{m}{\Sigma }}{b}_{\mathrm{cj}}=1.$

本发明具有的有益效果体现在：

数控装备类似于人体，其机能随着时间延长而逐渐衰退，具体表现为数控装备的功能与性能逐步劣化，导致数控装备发生无计划故障停机率增大，降低了数控装备运行可靠性。本发明专利通过对数控装备的性能劣化及可靠性变动进行监测、分析与预测，能够帮助企业发现数控装备潜在的失效模式，制定合理的维修计划以避免数控装备发生突发性失效，从而实现数控装备的可靠性增长，这对于提高数控装备利用率、减少数控装备维修费用、延长数控装备使用寿命以及构建数控装备健康状态监控体系等都具有重要的意义。

附图说明

图1本发明的总体原理图。

图2本发明的矢量量化图。

图3本发明的状态变迁拓扑结构图。

图4本发明的可靠性变动计算结果图。

具体实施方式

以下结合附图和实施例对本发明所述的技术方案作进一步介绍。

参照图1，一种基于离散隐马尔科夫链模型的数控装备可靠性变动计算方法主要包括数控装备的性能特征参数监测与提取、性能劣化模型、性能特征参数矢量量化、状态变迁模型和可靠性计算等模块。

1、性能特征参数监测与提取。

性能特征参数监测与提取是通过安装在数控装备上或周围的传感器等测试仪器来检测装备运行时的各种动态性能信号(如振动信号、运动信号、液压信号等)，经信号处理与分析提取出特征参数。所提取的性能特征参数与信号的类型及其处理方法相关，但是总体上它应遵循以下原则：(1)性能特征参数应具有明确的物理意义，(2)性能特征参数易于测量且相对稳定，(3)性能特征参数为反映数控装备可靠性变动的关键性参数等。例如，若用加速度传感器来采集数控装备的振动信号，并采用时域法来分析振动信号的时间历程，则所能够提取的性能特征参数有最大峰值、均方根指、峰值因子和裕度因子等。

根据事先布置的测试仪器，在数控装备使用过程中的l个给定时间点分别进行测量，得到动态性能信号。对采集的信号进行处理与分析，提取出能及时准确地反映数控装备性能状态变化的性能特征参数，即D＝{(T，Y)|(t_i，y_ij)；i＝1，2，…，l；j＝1，2，…，m}，其中t_i为第i个测量时间点，y_ij为在t_i时刻测量的第j个性能特征参数。

2、建立数控装备的性能劣化模型。

检验性能特征参数的离散时间序列{(t_i，y_ij)；i＝1，2，…，l；j＝1，2，…，m}的平稳性。若为非平稳，则需要采用d阶差分法将非平稳的{y_ij}时间序列转化为平稳的{w_ij}时间序列，即$w_{\mathrm{ij}}={▿}^d{y}_{\mathrm{ij}}.$

采用nARIMA算法对{w_ij}时间序列进行建模。由于第j个性能特征参数w_ij的取值不仅与其前p步的各个取值w_(i-1)j，w_(i-2)j，...，w_(i-p)j有关，而且与其前q步的随机扰动值ε_(i-1)j，ε_(i-2)j，…，ε_(i-q)j，则得到p阶自回归q阶滑动平均模型：

$w_{\mathrm{ij}}=\underset{z=1}{\overset{p}{\Sigma }}{\theta }_z{w}_{(i-z)j}-\underset{z=1}{\overset{q}{\Sigma }}{\phi }_z{ϵ}_{(i-z)j}+{ϵ}_{\mathrm{ij}}---\left(1\right)$

其中，θ_z为自回归系数，φ_z为滑动平均系数，${ϵ}_{\mathrm{ij}}~N(0,{\sigma }_{ϵ}^2)$为测量误差，为方差。引入后移算子B，则式(1)转换为式(2)：

$\left(\begin{array}{c}\theta \left(B\right)w_{\mathrm{ij}}=\phi \left(B\right){ϵ}_{\mathrm{ij}}\\ \theta \left(B\right)=1-\underset{z=1}{\overset{p}{\Sigma }}{\theta }_z{B}^z,\phi \left(B\right)=1-\underset{z=1}{\overset{q}{\Sigma }}{\phi }_z{B}^z\end{array}\right)---\left(2\right)$

其中，模型的阶次p和q以及未知参数θ₁，...，θ_p，φ₁，...，φ_q和σ_ε直接影响到模型精度。采用最小信息准则(AIC准则)来确定模型阶次，采用非线性最小二乘法来估计模型未知参数。一旦确定模型阶次和模型参数，则得到了非平稳nARIMA(p，d，q)性能劣化模型：

$y_{\mathrm{tj}}=\underset{z=1}{\overset{p}{\Sigma }}{\phi }_z{y}_{(t-z)j}+{ϵ}_{\mathrm{tj}}-\underset{z=1}{\overset{q}{\Sigma }}{\phi }_z{ϵ}_{(t-z)j}---\left(3\right)$

其中，φ_z为滑动平均系数，测量误差${ϵ}_{\mathrm{tj}}~N(0,{\sigma }_{ϵ}^2)$p为自回归阶次，表示模型预测值与前p阶性能参数值有关，q为滑动平均阶次，表示对模型预测值误差进行q阶的修正。

3、性能特征参数矢量量化

一旦建立了性能劣化模型，则能够预测数控装备在未来一段时间内的性能特征参数变化情况。例如，{y_tj，j＝1，…，m}表示由性能劣化模型预测所得的一组t时刻性能特征参数。这些性能特征参数从侧面反映了数控装备的隐含状态变化特征。

然而，由于数控装备本身非常复杂，其使用过程中所观测到的性能特征参数并不能与状态简单的一一对应。因此，本发明运用离散隐马尔科夫链模型，建立数控装备状态变迁模型来分析数控装备运行时的性能状态变化特征及其可靠性变动情况。

由于状态变迁模型的输入量必须为有限的离散值，而所获得的性能特征参数是一个连续实数值，为此采用Lloyd(劳埃德)算法对性能特征参数进行矢量量化处理以形成离散编码集合。

参照图2，性能特征参数的矢量量化过程是：根据分区向量(partition)和码本向量(codebook)，矢量量化将特征参数取值范围顺序划分成N个区域，每个区域映射一个离散观测值。矢量化所需预测时间点数T₀一般取三到六个单位时间点，在这T₀个预测时间点段中的第一到第T₀时间点，第j个特征参数值在每个区域的索引值定义为：

$v_1^j=1,v_N^j=N,$

$v_{i_0}^j=\mathrm{Index}\left(y_{t_{1}j}\right)=\left(\begin{array}{c}1,y_{t_{1}j}\le \mathrm{partition}(1,j)\\ i_0,\mathrm{partition}(i_0-1,j)\le y_{t_{1}j}\le \mathrm{partition}(i_0,j)\\ N,\mathrm{partition}(N-1,j)\le y_{t_{1}j}\end{array}\right)$

i₀＝2，…，N-1，t₁∈(1，…，T₀)；

其中，分区向量(partition)由装备的性能失效标准矩阵H来定义，即：

$H={\left[h_{i_{1}j}\right]}_{(N-1)\times m}=\left(\begin{array}{ccc}{h}_{11}& ···& h_{1m}\\ h_{21}& ···& h_{2m}\\ ···& ···& ···\\ h_{(N-1)1}& ···& h_{(N-1)m}\end{array}\right),$$\mathrm{partition}(i_1,j)=h_{i_{1}j},$表示初始定义的用于衡量第j个特征参数处于第i₁区域的界定值，i₁＝1，…，N-1

对预测所得的数控装备在1～T₀时刻的性能特征参数集{y_tj}进行矢量化，得到性能特征参数的输出观测序列：

$O^j=\{o_1^j,···,o_{T_0}^j\},$$o_{\tilde{t}}^j\in (v_1^j,,v_N^j),$$\tilde{t}=1,···,T_0,j=1,···,m$

其中，矢量化所需预测时间点T₀一般取第三到第六个单位时间点。

4、状态变迁模型

参照图3，建立装备的状态变迁模型。运行状态变迁概率是指装备运行状态不变或转向更差运行状态的概率。

在状态变迁拓扑结构图中，状态空间S＝{0，1，2，…，k，k+1}是数控装备性能运行状态集。运行状态包括正常状态、劣化状态和失效状态。图3中圆圈表示隐状态，{0}表示正常状态，{1，2，…，k}表示k个劣化状态，{k+1}表示失效状态；有向弧称为变迁弧，表示状态的转移。在不考虑维修的情况下，数控装备性能劣化过程是不可逆的；每个状态只能向其右侧更高编号的状态转移，同时每个状态也可以向自身转移。因此，令数控装备在t时刻的状态变迁概率矩阵A表示为：

$A={\left[a_{\mathrm{cd}}\right]}_{(k+1)\times (k+1)}=\left(\begin{array}{ccccc}{a}_{00}& a_{01}& ···& a_{0k}& a_{0(k+1)}\\ 0& a_{11}& ···& a_{1k}& a_{1(k+1)}\\ ···& ···& ···& ···& ···\\ 0& 0& ···& a_{\mathrm{kk}}& a_{k(k+1)}\\ 0& 0& ···& 0& a_{(k+1)(k+1)}\end{array}\right)$

其中，a_cd为数控装备在t时刻处于状态{c}而在t+1时刻处于状态{d}的概率，0≤a_cd≤1，且$\underset{d=0}{\overset{k+1}{\Sigma }}{a}_{\mathrm{cd}}=1.$

数控装备的状态由性能特征参数观测值来感知。根据数控装备在T0时刻的输出观测序列$O^j=\{o_1^j,···,o_{T_0}^j\},$令数控装备在t时刻的观测值概率矩阵B表示为：

$B={\left[b_{\mathrm{cj}}\right]}_{(k+1)\times m}=\left(\begin{array}{cccc}{b}_{01}& b_{02}& ···& b_{0m}\\ b_{11}& b_{12}& ···& b_{1m}\\ ···& ···& ···& ···\\ b_{k1}& b_{k2}& ···& b_{\mathrm{km}}\\ b_{(k+1)1}& b_{(k+1)2}& ···& b_{(k+1)m}\end{array}\right)$

其中，b_cj为数控装备在t时刻处于状态{c}时出现观测值$\left\{v_c^j\right\},c=0,···,k+1$的概率，0≤b_cj≤1，且$\underset{j=1}{\overset{m}{\Sigma }}{b}_{\mathrm{cj}}=1.$

初始条件下，数控装备处于正常运行状态，则初始状态的概率分布向量为π₀＝(1，0，…，0)，状态变迁概率矩阵A₀和观测值概率矩阵B₀采取随机方法选取，从而得到初始条件下数控装备的状态变迁模型λ₀＝(A₀，B₀，π。)。

对运行状态集、状态变迁概率矩阵A及观测值概率矩阵B定义后，采用离散隐马尔科夫链模型求解状态变迁概率矩阵A及观测值概率矩阵B，具体为：将性能特征参数在1～T₀时刻的输出观测序列$O^j=\{o_1^j,···,o_{T_0}^j\}$输入到初始状态变迁模型λ₀＝(A₀，B₀，π。)中，利用Baum-Welch(鲍姆-韦尔奇)算法进行迭代计算，使得模型参数逐步趋向更为合理的较优值，从而得到数控装备在t时刻的状态变迁模型λ＝(A，B，π)。为了验证该状态变迁模型的合理性，采用Forward-Backward(前后向)算法计算观测序列O在给定λ下的概率，即P(O|λ)。如果P(O|λ)超过期望值0.8，则认为得的状态变迁模型λ＝(A，B，π)是可行的。

5、可靠性计算

一旦确定数控装备在t时刻的状态变迁模型λ＝(A，B，π)，则计算数控装备的可靠性变动情况。具体过程是：令P_c(t)＝P(q_t＝c)表示数控装备在t时刻处于c状态的概率。根据切普曼-柯尔莫哥洛夫微分方程，有：

P′(t)＝P(t)·A   (4)

其中，P(t)＝(P₀(t)，P₁(t)，…，P_k(t)，P_(k+1)(t))为状态向量，P′(t)为P(t)的一阶微分状态向量，A为t时刻的状态变迁矩阵。对式(4)进行Laplace(拉普拉斯)变换，得：

$s·\left(\begin{array}{c}{P}_0\left(s\right)\\ P_1\left(s\right)\\ \text{·}\\ ·\\ ·\\ P_k\left(s\right)\\ P_{(k+1)}\left(s\right)\end{array}\right)-\left(\begin{array}{c}{P}_0\left(0\right)\\ P_1\left(0\right)\\ \text{·}\\ ·\\ ·\\ P_k\left(0\right)\\ P_{(k+1)}\left(0\right)\end{array}\right)\left(\begin{array}{ccccc}{a}_{00}& 0& ···& 0& 0\\ a_{01}& a_{11}& ···& 0& 0\\ ···& ···& ···& ···& ···\\ a_{0k}& a_{1k}& ···& a_{\mathrm{kk}}& 0\\ a_{0(k+1)}& a_{1(k+1)}& ···& a_{k(k+1)}& a_{(k+1)(k+1)}\end{array}\right)\left(\begin{array}{c}{P}_0\left(s\right)\\ P_1\left(s\right)\\ \text{·}\\ ·\\ ·\\ P_k\left(s\right)\\ P_{(k+1)}\left(s\right)\end{array}\right)---\left(5\right)$

若数控装备在初始条件通常处于正常状态，则有：

P(0)＝(P₀(0)，P₁(0)，…，P_k(0)，P_(k+1)(0))＝(1，0，…，0，0)

同时，将t时刻状态变迁矩阵A代入式(5)中，则计算得到：

P(s)＝(P₀(s)，P₁(s)，…，P_k(s)，P_(k+1)(s))

接着，对P(s)进行Laplace逆变换，则得到数控装备在t时刻处于不同状态的概率P(t)＝(P₀(t)，P₁(t)，…，P_k(t)，P_(k+1)(t))，计算出数控装备t时刻的可靠度R(t)＝1-P_(k+1)(t)等可靠性指标。进而，分析数控装备使用过程中的可靠性变动情况，发现数控装备的薄弱环节。

实施例：

该实施例给出了本发明在工程实践中的具体实施过程，同时验证了该发明的有效性。

参照图4，以欧泰OTM-650数控铣床为研究对象，通过模拟装备因反向间隙而引起的圆运动精度劣化，推断数控装备可靠性变动规律。试验装置包括欧泰OTM-650数控铣床、华中“世纪星”数控系统HMC-21M、海德汉KGM182平面光栅测量装置和普通PC机等。其中，平面光栅KGM182是一种用于非接触式测量装备二维运动轮廓精度的高精密测试仪器。通过NC编程使工作台在XOY平面内以半径为50mm和进给速度为2000mm/min作逆时针/顺时针圆运动，采用平面光栅在机测量装备的圆运动轨迹。让数控装备运行480个小时，每间隔大约20个小时采样1次，得到24组圆形轨迹。根据ISO230-4：1996给出的数控机床圆检验标准，检验与分析测量所得的24组圆形轨迹，计算出对应的圆滞后、圆偏差及半径偏差等圆运动精度特征量，如表1所示。

表1 圆运动精度特征量

由于每种特征量仅有24个实测值，为了提高拟合与预测精度，根据实测值采用插补再抽样法生成96个数据，共计有120个圆运动精度特征量。对于每种特征量，选取前100个数据作为模型训练集，后20个数据作为模型校验集。经检验发现圆运动精度特征量是非平稳的，为此采用一阶差分法将它们转化为平稳的时间序列{w_ij}。以圆滞后为例，根据AIC准则确定性能劣化模型的阶次，当(p，q)为(5，2)时，模型的信息函数达到最小值1.7939。进而，采用非线性最小二乘法来确定nARIMA(5，1，2)的模型参数，即：

$\left(\begin{array}{c}\theta \left(B\right)=1-{0.6847B}^1+{0.8273B}^2-{0.7267B}^3+{0.0877B}^4+{0.1747B}^2\\ \phi \left(B\right)=1-{0.0498B}^1-{0.9502B}^2\end{array}\right)$

从而，建立圆滞后劣化模型为：

y_ij＝0.6847y_(i-1)j-0.8273y_(i-2)j+0.7267y_(i-3)j-0.0877y_(i-4)j-0.1747y_(i-5)j+…

    +ε_ij-0.0498ε_(i-1)j-0.9502ε_(i-2)j；i＝5，…，100，j＝4

利用后20个数据来验证该模型，模型的FPE(最终预测偏差)为1.5459e-5。由此得出结论：该模型满足实际应用的要求。

类似建立圆滞后劣化模型的过程，可以分别建立如下半径偏差和圆偏差的性能劣化模型：

(1)半径偏差(F_max)劣化模型为：

y_ij＝y_(i-1)j+0.9165y_(i-2)j-1.43y_(i-3)j+0.5338y_(i-4)j+0.0033y_(i-5)j+…

     +ε_ij-0.0048ε_(i-1)j-0.9952ε_(i-2)j；i＝5，…，100，j＝1

(2)半径偏差(F_min)劣化模型为：

y_ij＝1.267y_(i-1)j-0.1327y_(i-2)j-0.723y_(i-3)j+0.4271y_(i-4)j-0.1187y_(i-5)j+0.2803y_(i-6)j+…

      +0.2942y_(i-7)j-0.1943y_(i-8)j+0.1601y_(i-9)j-0.1513y_(i-10)j+0.4735y_(i-11)j-0.166y_(i-12)j+…

      -0.4402y_(i-13)j+ε_ij-0.7509ε_(i-1)j+0.2458ε_(i-2)j；i＝13，…，100，j＝2

(3)圆偏差(G)劣化模型为：

y_ij＝0.5102y_(i-1)j+0.5102y_(i-2)j+ε_ij-0.0043ε_(i-1)j+0.207ε_(i-2)j；i＝2，…，100，j＝3

根据上述性能劣化模型，结合数控装备的加工工艺要求和使用规程，将数控装备失效状态划分为正常、轻微劣化、严重劣化和失效四个状态，并确定圆运动精度失效标准，如表2所示。接着，采用Lloyd算法对数控装备的性能特征参数进行矢量量化，得到数控装备t时刻的圆运动精度输出观察序列0。

表2 圆运动精度失效标准

由于数控装备在初始条件下处于正常运行状态，则初始状态的概率分布向量π₀＝[1，0，0，0]^T，状态变迁概率矩阵A₀和观测值概率矩阵B₀分别采取随机方法选取，有：

$A_0=\left(\begin{array}{cccc}0.3853& 0.2959& 0.2063& 0.1125\\ 0& 0.2198& 0.4064& 0.3738\\ 0& 0& 0.4989& 0.5011\\ 0& 0& 0& 1\end{array}\right),$$B_0=\left(\begin{array}{cccc}0.2343& 0.1668& 0.1296& 0.4693\\ 0.2575& 0.3889& 0.1219& 0.2317\\ 0.1961& 0.2752& 0.0932& 0.4355\\ 0.2710& 0.2810& 0.2115& 0.2365\end{array}\right)$

从而建立初始状态下的状态变迁模型λ₀＝(A₀，B₀，π₀)。将圆运动精度输出观察序列O代入λ₀＝(A₀，B₀，π₀)中，利用Baum-Welch算法，通过迭代计算来训练状态变迁模型。设置总的训练迭代次数为100，随着迭代次数的增加，模型输出对数概率值逐步增加，经过10次迭代后模型进入收敛饱和状态，输出概率值几乎不再增大，得到t时刻的状态变迁模型λ＝(A，B，π)，即：

$\pi =\left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right),$$A=\left(\begin{array}{cccc}0.9443& 0.0455& 0.0102& 0\\ 0& 0.6142& 0.3857& 0.0001\\ 0& 0& 0.9665& 0.0335\\ 0& 0& 0& 1\end{array}\right),$$B=\left(\begin{array}{cccc}1& 0& 0& 0\\ 1& 0& 0& 0\\ 1& 0& 0& 0\\ 0.1021& 0.4290& 0.1546& 0.3143\end{array}\right)$

同时，采用Forward-Backward算法来计算出t时刻圆运动精度输出观察序列O在λ＝(A，B，π)下的概率值P(O|λ)为80.84％。由此得出结论：建立的状态变迁模型λ＝(A，B，π)能够满足实际应用要求。

最后，将λ＝(A，B，π)中的π和A代入到式5中，求解该式，并经Laplace逆变换得到t时刻数控装备的状态向量P(t)＝(P₀(t)，P₁(t)，P₂(t)，P₃(t))，其中：

P₀(t)＝exp(0.9443t)

$P_1\left(t\right)=0.2757\exp \left(0.7793t\right)\frac{\exp \left(0.1651t\right)-\exp (-0.1651t)}{2}$

P₂(t)＝13.77exp(0.9443t)-13.93exp(0.9397t)+0.1633exp(0.6142t)

P₃(t)＝-8.28exp(0.9443t)+7.74exp(0.9397t)-0.0141exp(0.6142t)+0.5556exp(t)

则数控装备在t时刻的可靠度为R(t)＝1-P₃(t)，参照图4。由图4可知，数控装备运行大约300个小时之后可靠度已低于0.8，数控装备处于严重劣化状态，出现故障的概率很高。