融合寿命信息与多性能退化信息的元动作单元可靠性评估

Reliability Assessment of Meta-action Unit Integrating Lifetime Information and Multi-performance Degradation Information

摘要

针对机床可靠性建模困难、可靠性数据匮乏等问题，从元动作单元出发，提出一种融合寿命信息与多性能退化信息的可靠性评估方法 .利用“功能—运动—动作”（functionmotion-action,FMA）分解法将机床分解至最小运动单元，即元动作单元，机床的可靠性水平通过准确评估关键元动作单元的可靠性来保证.建立随机效应下非线性Wiener过程的性能退化模型，刻画试验样本退化过程存在的个体差异性、共同属性以及退化过程的非线性，在此基础上针对退化数据不足的问题，利用贝叶斯方法建立融合失效寿命信息的可靠性模型.考虑多性能退化失效相关性，基于Copula函数建立多性能退化可靠性评估模型，利用MCMC-Gibbs抽样算法实现参数估计.以某型号数控滚齿机的蜗杆转动元动作单元为研究对象，通过实例分析以及验证对比表明所提可靠性评估方法的可行性和优越性.

Abstract

To address the challenges of machine tool reliability modeling and the lack of reliability data， a reliability assessment method integrating lifetime information with multi-performance degradation information is proposed， starting from the basic motion units. Using the function-motion-action （FMA） decomposition method， the machine tool is broken down into its smallest motion units， i.e.， the meta-action units. The reliability level of the machine tool is ensured by accurately assessing the reliability of the key elemental action units. A non-linear Wiener process model under random effects is developed to capture individual variations， shared characteristics， and non-linearity in degradation processes. To address limited degradation data， a reliability model incorporating failure life information is established using Bayesian methods. Finally， considering correlations among multiple performance degradation failures， a comprehensive reliability assessment model is formulated using the Copula function， with parameter estimation achieved via the MCMC-Gibbs sampling algorithm. Focusing on the worm gear rotation meta-action unit of a specific computerized numerical control gear hobbing machine， case studies and comparative analysis highlight the method’s feasibility and superiority.

Graphical abstract

关键词

可靠性 / 信息融合 / 随机过程 / 元动作单元 / 贝叶斯方法

数控机床具有高复杂化和高寿命化的特性，如何在合理的时间与经济成本下准确了解其可靠性水平已经成为研究热点. 基于性能退化数据开展机床相关可靠性评估的方法被广泛使用. Wiener过程可以描述退化过程不单调的情况，且具有良好的计算特性，常被用于机床相关的可靠性分析中^［1-2］. 然而，传统方法中忽略随机影响因素和非线性退化特征的做法，往往与实际情况存在偏差. 目前，针对随机效应影响下的非线性随机过程研究仍主要围绕Wiener过程展开. Feng等^［3］考虑了性能退化过程不确定性和测量误差，推导了非线性退化场合的剩余寿命分布. Tang等^［4］研究了基于具有随机效应的非线性Wiener过程的加速退化过程建模问题，研究表明有必要将随机效应和非线性纳入加速退化建模中. 王汉禹等^［5］提出了基于核函数-Wiener过程退化模型，使用Wiener过程刻画轧辊退化趋势的强随机性特征，引入核函数捕捉轧辊的非线性退化路径. 张忠文等^［6］针对不同伺服驱动单元之间存在的个体差异特点，引入随机效应并给出了考虑个体差异的可靠性模型建立方案. 可以看出，针对随机效应下非线性Wiener过程的研究已经比较成熟，但基于随机过程建立可靠性模型也需要大量样本数据，对于机床或其部件而言，其各项性能指标的衰退周期较长且测试成本高，导致所获得的退化数据依旧是小规模的子样本^［7］.

针对数据不足的情况，目前研究主要集中在对数据层和算法层面的改进，如采用改进Bootstrap法扩充原始样本，构建可靠性模型，但这些方法所构建模型的精度无法保障，且方法适用范围较为有限. 因此需要充分利用所有的可靠性信息，在小子样本下构建高置信度的可靠性模型. 信息融合技术在数控机床可靠性研究领域的应用已经逐渐开展. 彭卫文^［8］考虑数控机床的多层次信息，提出了基于贝叶斯方法的独立来源层次客观数据的融合方法、基于贝叶斯网络的层次覆盖客观数据的融合方法以及基于贝叶斯网络的层次主客观信息的融合方法. 杨文^［9］提出了基于多源信号融合的数控铣床运行可靠性分析方法，通过对多传感器的数据进行数据预处理，实现数据层面的融合，随后将其作为输入数据，提取特征，根据基于状态距离的运行可靠性计算方法建立铣床的运行可靠性评估模型. 孙博^［10］以前期产品的可靠性信息和专家对现有产品可靠性水平判断结果为基础，通过贝叶斯融合方法与可靠性仿真信息融合，再利用改进的重要性重采样算法得到了伺服刀架在多源信息融合下的可靠性评估结果. 可以看出，目前关于机床多源信息融合下的可靠性评估大多都基于整机或部件层面的可靠性信息，但整机或部件所存在的固有问题，使能够获得的可靠性数据极为有限，以整机或部件为分析对象难以从数据中提取有效的退化信息，增加数据的不准确性，导致最终的评估结果出现较大偏差.

数控机床的功能和性能是通过其组成部件（零件）之间的相对运动予以实现和保障的^［11］. 基于机械结构运动和动力的传递，采用“功能-运动-动作”（FMA）的结构分解法将整机分解至最小运动单元（元动作单元）^［12-13］. Li等^［14］提出了基于元动作单元的维修决策. 张根保等^［15］提出一种基于非线性Gamma随机过程和混合 Copula 函数的元动作单元性能可靠性分析模型. 鞠萍华等^［16］从元动作单元出发，建立了关于元动作单元的故障概率模型. 元动作单元是机床中最小的结构单元，其可靠性分析是整机可靠性的基础，以关键元动作单元为研究对象，可以避免结构分析困难，试验难以进行，数据难以收集，成本高，周期长等问题.

另外，元动作单元运行过程中，故障是在多个因素共同作用下发生的，针对具有多个性能退化指标的元动作单元，其退化失效之间具有相关性. Zhu等^［17］基于电池退化数据构建了二元Wiener过程退化模型. 董庆来等^［18］考虑两类竞争失效形式，构建了基于二元Wiener过程的可靠性模型. 多元退化过程通常依赖于线性假设，并且退化模型存在结构复杂、求解困难的问题. 而Copula函数能够对非线性依赖进行捕捉建模，且更容易估计和拟合，这在处理复杂的实际问题时尤其重要.

综上所述，本文提出一种融合寿命信息与多性能退化信息的可靠性评估方法. 以元动作单元为研究对象，首先用考虑个体差异和共同属性的非线性Wiener过程描述其性能退化过程，建立性能退化模型；其次将历史失效寿命数据作为先验信息，利用贝叶斯方法融合性能退化信息和寿命信息，构建融合了失效寿命信息的单性能退化可靠性模型；然后根据赤池信息准则（Akaike information criterion，AIC）和贝叶斯信息准则（Bayesian information criterion，BIC）选取合适的Copula函数，用Copula函数联合多个边缘失效分布，得到考虑多退化失效相关性的可靠性评估模型，模型中的未知参数通过MCMC-Gibbs抽样法逐步获得；最后以某型号数控滚齿机中的蜗杆转动元动作单元为例，验证本文方法的可行性与优越性.

1 理论介绍

1.1 FMA分解

根据元动作理论，采用FMA结构化分解方法，将机械产品的功能一步步分解到运动的最小运动单元——元动作. FMA结构化分解示意图如图1所示，图中MAC为元动作链（meta-action chain）.

1.2 元动作单元

元动作单元（meta-action unit， MAU）是实现元动作的基本结构单元. 元动作单元比部件更细化，可用信息更多，且它可以自成体系，体现出零部件之间的机械作用，具有更强的综合性，元动作单元在结构上独立，可以独立进行设计、分析和试验^［12］. MAU的结构组成如图2所示.

相比于整机及其部件，对元动作单元进行可靠性试验的可行性更高，元动作单元的可靠性信息更容易获得，且由于其结构清晰，运动形式单一，从收集到的单元数据中可以获取更直观有效的性能退化信息，为可靠性评估工作奠定了良好的基础.

2 基于贝叶斯的信息融合可靠性模型

2.1 随机效应下非线性Wiener过程的可靠性建模

2.1.1 随机效应下非线性Wiener过程的性能退化模型

由于正态分布方差的先验分布一般被认为是逆伽马分布，则假设{\sigma }^{\mathrm{2}}~IGamma（a，b），在给定{\sigma }^{\mathrm{2}}的条件下，参数μ服从正态分布μ|{\sigma }^{\mathrm{2}}~N（{\theta }_{\mu },{\lambda }_{\mu }/{\sigma }^{\mathrm{2}}）. 根据文献［19］，考虑产品个体差异性和共同属性的非线性Wiener过程模型可表示为：

\left\{\begin{array}{l}Y(t)=Y(\mathrm{0})+\mu \Lambda (t)+\sigma B\left[\Lambda (t)\right],t\ge \mathrm{0}\\ {\sigma }^{\mathrm{2}}\sim \mathrm{I}\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}(a,b)\\ \mu \left|{\sigma }^{\mathrm{2}}\sim N({\theta }_{\mu },{\lambda }_{\mu }/{\sigma }^{\mathrm{2}})\right.\end{array}\right.（1）

式中：μ表示退化速率的漂移系数；σ表示扩散系数；B\left[\Lambda (t)\right]为布朗运动，用来描述性能退化量的不确定性. 目前针对描述非线性特征的连续非线性函数\Lambda (t)有两种模型：{\Lambda }_{\mathrm{1}}(t)={\Lambda }_{\mathrm{1}}(t,d)=t^d；{\Lambda }_{\mathrm{2}}(t)={\Lambda }_{\mathrm{2}}(t,d)=\mathrm{e}\mathrm{x}\mathrm{p}(dt)，应当根据退化数据选用合适的转换函数.

Wiener过程退化值变化量的分布为\mathrm{\Delta }{y}_{ij}\left|\right.\mu ,{\sigma }^{\mathrm{2}}\sim N(\mu \mathrm{\Delta }{\Lambda }_i(t_j),{\sigma }^{\mathrm{2}}\mathrm{\Delta }{\Lambda }_i(t_j))，Y_{ij}表示第i个产品在t_j时刻的退化测量值(i=\mathrm{1,2},\cdots ,N;j=\mathrm{1,2},\cdots ,M)，\mathrm{\Delta }{y}_{ij}=Y_{{}_i}(t_j)-Y_{{}_i}(t_{j-\mathrm{1}})，\mathrm{\Delta }{\Lambda }_i(t_j)={\Lambda }_i(t_j)-{\Lambda }_i(t_{j-\mathrm{1}})，它的概率密度函数为：

f(\mathrm{\Delta }{y}_{ij}\left|\mu ,{\sigma }^{\mathrm{2}}\right.)=\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{{\left[\mathrm{\Delta }{y}_{ij}-\mu \mathrm{\Delta }\Lambda (t)\right]}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}\mathrm{\Delta }\Lambda (t)}\right\}}{\sqrt[]{\mathrm{2}\mathrm{\pi }{\sigma }^{\mathrm{2}}\mathrm{\Delta }\Lambda (t)}}（2）

根据全概率公式，考虑产品个体差异性和共同属性下非线性Wiener过程的概率密度函数，得

f_{Y(t)}(\mathrm{\Delta }{y}_{ij})=\frac{\mathrm{\Gamma }\left(a+\frac{\mathrm{1}}{\mathrm{2}}\right){\left[\mathrm{1}+\frac{{\left(\mathrm{\Delta }y-{\theta }_{\mu }\mathrm{\Delta }\Lambda \right)}^{\mathrm{2}}}{\mathrm{2}b\left({\lambda }_{\mu }\mathrm{\Delta }{\Lambda }^{\mathrm{2}}+\mathrm{\Delta }\Lambda \right)}\right]}^{-a-\mathrm{1}/\mathrm{2}}}{\sqrt[]{\mathrm{2}\mathrm{\pi }b\left({\lambda }_{\mu }\mathrm{\Delta }{\Lambda }^{\mathrm{2}}+\mathrm{\Delta }\Lambda \right)}\mathrm{\Gamma }(a)}（3）

式中：\mathrm{\Gamma }(a)为Gamma函数. 基于随机效应的非线性Wiener过程的似然函数为：

L={\displaystyle \prod _{i=\mathrm{1}}^N}{\displaystyle \prod _{j=\mathrm{1}}^M}\frac{\mathrm{\Gamma }\left(a+\frac{\mathrm{1}}{\mathrm{2}}\right){\left[\mathrm{1}+\frac{{\left(\mathrm{\Delta }y-{\theta }_{\mu }\mathrm{\Delta }\Lambda \right)}^{\mathrm{2}}}{\mathrm{2}b\left({\lambda }_{\mu }\mathrm{\Delta }{\Lambda }^{\mathrm{2}}+\mathrm{\Delta }\Lambda \right)}\right]}^{-a-\mathrm{1}/\mathrm{2}}}{\sqrt[]{\mathrm{2}\mathrm{\pi }b\left({\lambda }_{\mu }\mathrm{\Delta }{\Lambda }^{\mathrm{2}}+\mathrm{\Delta }\Lambda \right)}\mathrm{\Gamma }(a)}（4）

可以看出该似然函数的形式十分复杂，难以通过极大似然法得到退化模型的未知参数的解析解.

2.1.2 随机效应下非线性Wiener过程的可靠性模型

令C表示性能退化过程的失效阈值，产品失效时间T为性能退化量首次达到失效阈值C的时间，则T表示为：

T=\mathrm{i}\mathrm{n}\mathrm{f}\{t\left|Y(t)\ge C,t>\mathrm{0}\}\right.（5）

已知一般线性Wiener过程所对应的首次达到失效阈值时间（后续简称首达时间）分布为逆高斯分布，其概率密度函数为：

f(t)=\frac{C}{\sqrt[]{\mathrm{2}\mathrm{\pi }{\sigma }^{\mathrm{2}}{t}^{\mathrm{3}}}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{(C-\mu t)}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}t}\right]（6）

若退化过程不存在个体差异，即参数\{\mu ,{\sigma }^{\mathrm{2}}\}是固定参数，但是存在非线性特性，根据式（6）可以推出固定效应下的非线性Wiener过程所对应的失效概率密度函数为：

f(t)=\frac{C_{\mathrm{W}}}{\sqrt[]{\mathrm{2}\mathrm{\pi }\Lambda {(t)}^{\mathrm{3}}}\sigma }\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{{\left[C_{\mathrm{W}}-\mu \Lambda (t)\right]}^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}\Lambda (t)}\right\}\Lambda \mathrm{\text{'}}(t)（7）

当样本之间考虑个体差异和共同属性时，参数\{\mu ,{\sigma }^{\mathrm{2}}\}为随机变量，则根据随机变量的全概率公式可以得到产品的失效概率密度为：

f_{\mathrm{T}}(t)=\frac{\mathrm{\Gamma }\left(a+\frac{\mathrm{1}}{\mathrm{2}}\right)C_{\mathrm{W}}{\left\{\mathrm{1}+\frac{{\left[C_{\mathrm{W}}-{\theta }_{\mu }\Lambda (t)\right]}^{\mathrm{2}}}{\mathrm{2}b\left[{\lambda }_{\mu }\Lambda {(t)}^{\mathrm{2}}+\Lambda (t)\right]}\right\}}^{-a-\mathrm{1}/\mathrm{2}}}{\sqrt[]{\mathrm{2}\mathrm{\pi }\Lambda {(t)}^{\mathrm{3}}\left\{b\left[({\lambda }_{\mu }\Lambda (t)+\mathrm{1}\right]\right\}\mathrm{\Gamma }(a)}}（8）

其失效分布函数为：

F_{\mathrm{T}}(t)={\int }_{\mathrm{0}}^t{f}_{\mathrm{T}}(t)\mathrm{d}t=T_{\mathrm{2}a}\left[\sqrt[]{\frac{a}{b}}\frac{{\theta }_{\mu }\Lambda (t)-C_{\mathrm{W}}}{\sqrt[]{{\lambda }_{\mu }\Lambda {(t)}^{\mathrm{2}}+\Lambda (t)}}\right]（9）

式中：T_{\mathrm{2}a}为自由度为\mathrm{2}a的T分布函数.

根据可靠度函数与失效分布函数之间的关系可得，随机效应下非线性Wiener过程的可靠度函数为：

\begin{array}{l}R(t\left|{\theta }_{\mu },{\lambda }_{\mu },a,b,d)=\mathrm{1}-\right.\\ T_{\mathrm{2}a}\left[\sqrt[]{\frac{a}{b}}\frac{{\theta }_{\mu }\Lambda (t)-C_{\mathrm{W}}}{\sqrt[]{{\lambda }_{\mu }\Lambda {(t)}^{\mathrm{2}}+\Lambda (t)}}\right]\end{array}（10）

2.2 基于失效寿命数据的先验分布确定

由2.1.1节可知，性能退化模型的未知参数集为\Psi =\left\{{\theta }_{\mu },{\lambda }_{\mu },a,b,d\right\}. 文献［20］研究了随机效应下隐含非线性Wiener退化过程模型的参数估计的性质，基于无偏估计分析了未知参数估计的精度与其影响因素的关系. 经分析可得，未知参数集\Psi中，\left\{a\right.,b,d\}这3个参数的估计精度与进行模型参数估计时使用的数据量成正比，数据量的大小是估计精度的决定性因素；而未知参数\{{\theta }_{\mu },{\lambda }_{\mu }\}的估计精度除了与数据量有正比关系之外还会受到观测时长的影响，观测时间越长，精度越高，因此在无法得到有效的先验性能退化信息或缺失信息的情况下，可以通过融合失效寿命数据来提高未知参数的估计精度. 综合考虑各个未知参数估计的性质以及贝叶斯方法的计算复杂程度，可以将未知参数\left\{a\right.,b,d\}视为固定参数，它们的先验分布为无信息先验分布，未知参数\{{\theta }_{\mu },{\lambda }_{\mu }\}的先验分布将根据失效寿命数据和最大熵法求得. 确定随机效应下非线性Wiener过程模型中未知参数\{{\theta }_{\mu },{\lambda }_{\mu }\}的先验分布过程如下.

假设固定参数\left\{a\right.,b,d\}已给定，令失效寿命数据T_{\mathrm{W}}=\{T_{\mathrm{W};\mathrm{1}},T_{\mathrm{W};\mathrm{2}},\cdots ,T_{\mathrm{W};N}\}，T_{\mathrm{W};N\mathrm{\text{'}}}表示第N\mathrm{\text{'}}个产品在退化过程中首次达到失效阈值C的时间，N\mathrm{\text{'}}\in [\mathrm{1},N]. 假设\Lambda (t)=t^d，由非线性Wiener过程的性质可知^［20］，失效时间T_{\mathrm{W};N\mathrm{\text{'}}}服从逆高斯分布，可以推导出T_{\mathrm{W};N\mathrm{\text{'}}}的似然函数为：

\begin{array}{l}L\left({\mu }_{N\mathrm{\text{'}}}\left|T_{\mathrm{W};N\mathrm{\text{'}}}\right.\right)=\frac{C_{\mathrm{W}}-{\mu }_{N\mathrm{\text{'}}}\left(T_{{}_{\mathrm{W};N\mathrm{\text{'}}}}^{{}^{\widehat{d}}}-\widehat{d}{T}_{{}_{\mathrm{W};N\mathrm{\text{'}}}}^{{}^{\widehat{d}}}\right)}{\sqrt[]{\mathrm{2}\mathrm{\pi }{T}_{{}_{\mathrm{W};N\mathrm{\text{'}}}}^{{}^{\mathrm{3}}}}\widehat{\sigma }}\times \\ \mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\left(C_{\mathrm{W}}-{\mu }_{N\mathrm{\text{'}}}{T}_{{}_{\mathrm{W};N\mathrm{\text{'}}}}^{{}^{\widehat{d}}}\right)}^{\mathrm{2}}}{\mathrm{2}{T}_{{}_{\mathrm{W};N\mathrm{\text{'}}}}^{{}^{\widehat{d}}}{\widehat{\sigma }}^{\mathrm{2}}}\right]\end{array}（11）

式中：{\mu }_{N\mathrm{\text{'}}}为第N\mathrm{\text{'}}个产品首达时间分布所对应的漂移参数.

在给定失效寿命数据T_{\mathrm{W}}=\left\{T_{W;\mathrm{1}},T_{W;\mathrm{2}},\right.\cdots ,\left.T_{\mathrm{W};N}\right\}的条件下，利用MATLAB中的fminsearch函数求取参数{\mu }_{N\mathrm{\text{'}}}的极大似然估计值，得到的结果表示为{\widehat{\mu }}_{N\mathrm{\text{'}}}=\{{\widehat{\mu }}_{N{\mathrm{\text{'}}}_{\mathrm{1}}},{\widehat{\mu }}_{N{\mathrm{\text{'}}}_{\mathrm{2}}},\cdots ,{\widehat{\mu }}_{N{\mathrm{\text{'}}}_N}\}，由此便可以推出未知参数\{{\theta }_{\mu },{\lambda }_{\mu }\}的先验分布均值，如式（12）所示.

\left\{\begin{array}{l}{\overline{\theta }}_{\mu }=\frac{\mathrm{1}}{N}{\displaystyle \sum _{N\mathrm{\text{'}}=\mathrm{1}}^N}{\widehat{\mu }}_{N\mathrm{\text{'}}}\\ {\overline{\lambda }}_{\mu }=\sqrt[]{\frac{\mathrm{1}}{N-\mathrm{1}}{\displaystyle \sum _{N\mathrm{\text{'}}=\mathrm{1}}^N}({\widehat{\mu }}_{N\mathrm{\text{'}}}-{\overline{\theta }}_{\mu }{)}^{\mathrm{2}}}\end{array}\right.（12）

在求得了未知参数\{{\theta }_{\mu },{\lambda }_{\mu }\}的先验分布的均值之后，需要找合适的分布，使得该分布的均值为此值并且包含最少的先验分布信息，考虑到只有未知参数的先验均值这一项已知信息，采用最大熵法构建先验分布最为合适. 在给定约束条件下，不确定性最大的分布即为该随机变量的最大熵分布^［21］. 根据已有信息，目前只得到了未知参数\{{\theta }_{\mu },{\lambda }_{\mu }\}的一阶矩（先验均值），因此参数{\theta }_{\mu }的熵函数为：

H[\pi ({\theta }_{\mu })]=-\int \pi ({\theta }_{\mu })\mathrm{l}\mathrm{n}[\pi ({\theta }_{\mu })]\mathrm{d}{\theta }_{\mu }（13）

\pi ({\theta }_{\mu })约束条件为：

\left\{\begin{array}{l}\int \pi ({\theta }_{\mu })\mathrm{d}{\theta }_{\mu }=\mathrm{1}\\ \int {\theta }_{\mu }\pi ({\theta }_{\mu })\mathrm{d}{\theta }_{\mu }={\overline{\theta }}_{\mu }\end{array}\right.（14）

要想联立公式求解，可以采用拉格朗日乘子法构造等式方程，求得原函数中各个变量的解. 构造函数：

F_L=H[\pi ({\theta }_{\mu })]+g_{\mathrm{1}}\left[\int \pi ({\theta }_{\mu })d{\theta }_{\mu }-\mathrm{1}\right]+

g_{\mathrm{2}}\left[\int {\theta }_{\mu }\pi ({\theta }_{\mu })\mathrm{d}{\theta }_{\mu }-{\overline{\theta }}_{\mu }\right]（15）

在偏导数\frac{\partial F_L}{\partial {\theta }_{\mu }}=\mathrm{0}的条件下求得极值点：

\pi ({\theta }_{\mu })=\mathrm{e}\mathrm{x}\mathrm{p}(g_{\mathrm{2}}{\theta }_{\mu }+g_{\mathrm{1}})（16）

可以解得：

\left\{\begin{array}{l}{g}_{\mathrm{1}}=-\mathrm{l}\mathrm{n}{\overline{\theta }}_{\mu }\\ g_{\mathrm{2}}=-\mathrm{1}/{\overline{\theta }}_{\mu }\end{array}\right.（17）

至此，依据最大熵原理可以确定参数{\theta }_{\mu }的先验分布为：

\pi ({\theta }_{\mu })=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\theta }_{\mu }}{{\overline{\theta }}_{\mu }}-\mathrm{l}\mathrm{n}{\overline{\theta }}_{\mu }\right)（18）

同理可得，未知参数{\lambda }_{\mu }的先验分布为：

\pi ({\lambda }_{\mu })=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\lambda }_{\mu }}{{\overline{\lambda }}_{\mu }}-\mathrm{l}\mathrm{n}{\overline{\lambda }}_{\mu }\right)（19）

2.3 寿命信息与性能退化信息融合可靠性模型

结合观测数据，对模型中的未知参数进一步更新，即确定未知参数的后验分布.

令模型的联合先验分布为\pi ({\theta }_{\mu },{\lambda }_{\mu })，认为参数{\theta }_{\mu }和{\lambda }_{\mu }之间相互独立，则联合先验分布为\pi ({\theta }_{\mu },{\lambda }_{\mu })=\pi ({\theta }_{\mu })\cdot \pi ({\lambda }_{\mu }). 根据贝叶斯理论，模型未知参数的联合后验分布可以表示为：

\pi ({\theta }_{\mu },{\lambda }_{\mu }\left|\right.Y)\propto \pi ({\theta }_{\mu },{\lambda }_{\mu })L(Y\left|{\theta }_{\mu },{\lambda }_{\mu },a,b,d\right.)=

\pi ({\theta }_{\mu },{\lambda }_{\mu }){\displaystyle \prod _{i=\mathrm{1}}^N}{\displaystyle \prod _{j=\mathrm{1}}^M}\frac{\mathrm{\Gamma }\left(a+\frac{\mathrm{1}}{\mathrm{2}}\right){\left[\mathrm{1}+\frac{{\left(\mathrm{\Delta }y-{\theta }_{\mu }\mathrm{\Delta }\Lambda \right)}^{\mathrm{2}}}{\mathrm{2}b\left({\lambda }_{\mu }\mathrm{\Delta }{\Lambda }^{\mathrm{2}}+\mathrm{\Delta }\Lambda \right)}\right]}^{-a-\mathrm{1}/\mathrm{2}}}{\sqrt[]{\mathrm{2}\pi b\left({\lambda }_{\mu }\mathrm{\Delta }{\Lambda }^{\mathrm{2}}+\mathrm{\Delta }\Lambda \right)}\mathrm{\Gamma }(a)}（20）

求得待估参数的后验分布\pi ({\theta }_{\mu },{\lambda }_{\mu }\left|\right.Y)之后，便可以得到融合单性能退化信息与寿命信息的可靠度函数为：

R(t)={\int }_{\mathrm{0}}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\pi ({\theta }_{\mu },{\lambda }_{\mu }\left|Y\right.)R(t\left|{\theta }_{\mu },{\lambda }_{\mu },a,b,d)\mathrm{d}{\theta }_{\mu }\mathrm{d}{\lambda }_{\mu }\right.（21）

同时也可以推出第i个产品在未来时间节点t_{i,M_i+\mathrm{1}}的性能退化的预测值为：

f\left[Y(t_{i,M_i+\mathrm{1}})\left|Y\right.\right]={\int }_{\mathrm{0}}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\pi ({\theta }_{\mu },{\lambda }_{\mu }\left|Y\right.)\times

f_{Y(t)}(\mathrm{\Delta }{y}_{ij}\left|{\theta }_{\mu },{\lambda }_{\mu },a,b,d\right.)\mathrm{d}{\theta }_{\mu }\mathrm{d}{\lambda }_{\mu }（22）

基于贝叶斯的寿命信息与性能退化信息融合可靠性模型的建模流程如图3所示.

3 基于Copula函数的多性能退化可靠性模型

假设联合分布函数H具有多个边界分布F_{\mathrm{1}}(x_{\mathrm{1}})，F_{\mathrm{2}}(x_{\mathrm{2}})，…，F_L(x_L)，则存在Copula函数C，对任意的x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_L，x_L\in R，有

H\left(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_L\right)=C\left[F_{\mathrm{1}}\left(x_{\mathrm{1}}\right),F_{\mathrm{2}}\left(x_{\mathrm{2}}\right),\cdots ,F_L\left(x_L\right)\right]（23）

联合概率密度函数h\left(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_L\right)为^［22］：

h\left(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_L\right)=\frac{{\partial }^{L}C\left(u_{\mathrm{1}},u_{\mathrm{2}},\cdots ,u_L\right)}{\partial u_{\mathrm{1}},\partial u_{\mathrm{2}},\cdots ,\partial u_L}{\displaystyle \prod _{l=\mathrm{1}}^L}\frac{\partial F_l\left(x_l\right)}{\partial x_l}=

c\left(u_{\mathrm{1}},u_{\mathrm{2}},\cdots ,u_L\right){\displaystyle \prod _{l=\mathrm{1}}^L}{f}_l\left(x_l\right)（24）

基于构造方法的不同Copula函数有很多种，常用二元Copula函数的性质及其秩相关系数计算式见参考文献［23］. 需要通过模型评价准则来进行Copula函数选取. 赤池信息准则（AIC）和贝叶斯信息准则（BIC）是最常用的Copula函数评价准则，AIC和BIC的计算值A_{\mathrm{I}\mathrm{C}}、B_{\mathrm{I}\mathrm{C}}分别为：

A_{\mathrm{I}\mathrm{C}}=-\mathrm{2}{L}_{\mathrm{C}}(\widehat{\theta })+\mathrm{2}{k}_{\mathrm{c}}=-\mathrm{2}{\displaystyle \sum _{i_{\mathrm{c}}=\mathrm{1}}^{m_{\mathrm{c}}}}\mathrm{l}\mathrm{n}c({\mu }_{{\mathrm{1}}_{i_{\mathrm{c}}}},{\mu }_{{\mathrm{2}}_{i_{\mathrm{c}}}};\widehat{\theta })+\mathrm{2}{k}_{\mathrm{c}}（25）

B_{\mathrm{I}\mathrm{C}}=-\mathrm{2}{L}_{\mathrm{C}}(\widehat{\theta })+k_{\mathrm{c}}\mathrm{l}\mathrm{n}{m}_{\mathrm{c}}=

-\mathrm{2}{\displaystyle \sum _{i_{\mathrm{c}}=\mathrm{1}}^{m_{\mathrm{c}}}}\mathrm{l}\mathrm{n}c({\mu }_{{\mathrm{1}}_{i_{\mathrm{c}}}},{\mu }_{{\mathrm{2}}_{i_{\mathrm{c}}}};\widehat{\theta })+k_{\mathrm{c}}\mathrm{l}\mathrm{n}{m}_{\mathrm{c}}（26）

式中：L_{\mathrm{C}}(\widehat{\theta })为模型的极大似然函数值；k_{\mathrm{c}}为模型参数的数量；m_{\mathrm{c}}为样本数量.

为建立多性能退化模型，令\left\{Y_{\mathrm{1}}(t),Y_{\mathrm{2}}(t),\right.\cdots ,\left.Y_L(t)\right\}表示产品在t时刻的多个性能退化量，对应的失效阈值记为\{C_{\mathrm{1}},C_{\mathrm{2}},\cdots ,C_L\}，失效时间记为\{T_{\mathrm{1}},T_{\mathrm{2}},\cdots ,T_L\}，其各自的失效分布函数记为\{F_{T_{\mathrm{1}}}(t),F_{T_{\mathrm{2}}}(t),\cdots ,F_{T_L}(t)\}，产品寿命可以表示为：T=\mathrm{m}\mathrm{i}\mathrm{n}(T_{\mathrm{1}},T_{\mathrm{2}},\cdots ,T_L). 多性能退化过程的可靠度函数为：

R(t)=P(T>t)=P\left[\mathrm{m}\mathrm{i}\mathrm{n}\left(T_{\mathrm{1}},T_{\mathrm{2}},\cdots ,T_L\right)>t\right]=

P\left(T_{\mathrm{1}}>t,T_{\mathrm{2}}>t,\cdots ,T_L>t\right)（27）

若存在多个性能退化量相互关联的情况，可以用Copula函数来描述它们之间的相关关系，产品的联合分布函数为：

H(t_{\mathrm{1}},t_{\mathrm{2}},\cdots ,t_L)=C\left[{F_T}_{{}_{\mathrm{1}}}\left(t_{\mathrm{1}}\right),{F_T}_{{}_{\mathrm{2}}}\left(t_{\mathrm{2}}\right),\cdots ,{F_T}_{{}_L}\left(t_L\right);\theta \right]（28）

已知可靠度R(t)=\mathrm{1}-F(t)，则产品考虑多元性能退化失效相关性的可靠性模型为：

R^C(t)=P(T>t)=P(T_{\mathrm{1}}>t_{\mathrm{1}},T_{\mathrm{2}}>t_{\mathrm{2}},\cdots ,T_L>t_L)=

\mathrm{1}-{\displaystyle \sum _{i=\mathrm{1}}^L}{F}_{T_i}(t)+\sum _{\mathrm{1}\le i\le j\le L}C\left[{F_T}_{{}_i}\left(t_i\right),{F_T}_{{}_j}\left(t_j\right)\right]+\cdots +

{(-\mathrm{1})}^{L}C\left[{F_T}_{{}_{\mathrm{1}}}\left(t_{\mathrm{1}}\right),{F_T}_{{}_{\mathrm{2}}}\left(t_{\mathrm{2}}\right),\cdots ,{F_T}_{{}_L}\left(t_L\right)\right]（29）

以二元性能退化下失效相关的可靠性模型为例：

R^C(t)=R_{\mathrm{1}}(t)+R_{\mathrm{2}}(t)-\mathrm{1}+C\left[F_{T_{\mathrm{1}}}(t),F_{T_{\mathrm{2}}}(t);\theta \right]（30）

4 模型参数估计与可靠性评估

由上述建立的可靠性模型可以看出，随机效应下非线性Wiener过程的似然函数与贝叶斯融合模型的后验分布含有多个未知参数，并且函数的形式也十分复杂. 马尔可夫链蒙特卡罗（Markov chain Monte Carlo， MCMC）算法是一种基于统计模拟的参数估计方法，广泛应用于函数形式复杂，无法直接计算的情况. MCMC的核心思想是通过构建一个马尔可夫链，使其平稳分布与目标分布一致，通过足够长时间的模拟，马尔可夫链的状态分布将趋近目标分布，从而可以通过采样得到目标分布的近似样本，进而用于参数估计或其他统计推断. MCMC-Gibbs抽样是MCMC方法的一种特殊形式，特别适用于多参数情况下的联合分布抽样. 鉴于此，本文采用MCMC-Gibbs抽样方法进行估计. MCMC-Gibbs抽样的基本步骤参考文献［24］.

在使用MCMC-Gibbs抽样算法求得未知参数\left\{{\theta }_{\mu },{\lambda }_{\mu },a,b,d\right\}的估计值之后，对Copula函数的相关参数\theta进行估计，参数\theta可以采用极大似然估计（maximum likelihood estimation，MLE）法求取. 记各个性能退化量的边缘失效分布函数为F_{T_i}(t)，i=\mathrm{1,2},\cdots ,L，根据式（24），基于Copula函数的性能退化量的联合边缘失效概率密度函数为：

f_T=c\left[F_{T_{\mathrm{1}}}\left(t\right),F_{T_{\mathrm{2}}}\left(t\right),\cdots ,F_{T_L}\left(t\right);\theta \right]{\displaystyle \prod _{i=\mathrm{1}}^L}{f}_{T_i}(t)（31）

以二元性能退化为例，令{\Psi }_{\mathrm{1}}=\left\{{{\theta }_{\mu }}_{{}_{\mathrm{1}}},{{\lambda }_{\mu }}_{{}_{\mathrm{1}}},\right.a_{\mathrm{1}},\left.b_{\mathrm{1}},d_{\mathrm{1}}\right\}和{\Psi }_{\mathrm{2}}=\left\{{{\theta }_{\mu }}_{{}_{\mathrm{2}}},{{\lambda }_{\mu }}_{{}_{\mathrm{2}}},a_{\mathrm{2}},b_{\mathrm{2}},d_{\mathrm{2}}\right\}分别表示两种性能退化下失效分布模型的未知参数，其对数似然函数可以表示为：

\mathrm{l}\mathrm{n}\left[L(\theta \left|{\Psi }_{\mathrm{1}},{\Psi }_{\mathrm{2}})\right.\right]=

{\displaystyle \sum _{i=\mathrm{1}}^N}\left[\mathrm{l}\mathrm{n}\left(c\left(F_{T_{\mathrm{1}}}(t_i;{\Psi }_{\mathrm{1}}),F_{T_{\mathrm{2}}}(t_i;{\Psi }_{\mathrm{2}});\theta \right)\right)\right.+

\left.\mathrm{l}\mathrm{n}\left(f_{T_{\mathrm{1}}}(t_i);{\Psi }_{\mathrm{1}}\right)+\mathrm{l}\mathrm{n}\left(f_{T_{\mathrm{2}}}(t_i);{\Psi }_{\mathrm{2}}\right)\right]（32）

将已经求得的未知参数\left\{{\theta }_{\mu },{\lambda }_{\mu },a,b,d\right\}的估计值代入上述对数似然函数中，利用MATLAB中的fminsearch函数便可以求出Copula函数的参数\theta的估计值. 至此便可以得到全部的未知参数估计值. 融合寿命信息与多性能退化信息可靠性评估方法流程图如图4所示.

5 实 例

5.1 某型号数控滚齿机FMA分解

FMA分解按照“功能-运动-动作”对滚齿机进行分解，分解图如图5所示. 通过对滚齿机的可靠性分析得到该滚齿机中关键元动作单元为蜗杆转动元动作单元，考虑保密要求，不进行具体阐述.

5.2 蜗杆转动元动作单元可靠性评估

5.2.1 可靠性模型

搭建蜗轮蜗杆副可靠性试验台，周期性监测收集齿侧间隙和蜗杆角位移的离散退化数据，试验台主要包括试验台架、蜗轮蜗杆、蜗轮轴、传感器和磁粉制动器，试验台架上安装电机支架，电机支架上固定电机，电机通过花键轴连接蜗杆，蜗杆蜗轮啮合.蜗轮蜗杆副可靠性试验台如图6所示.

结合失效寿命数据，使用上述方法完成可靠性建模，得到蜗杆转动元动作单元的可靠性评估结果. 如图7和图8分别为在蜗杆转动元动作单元可靠性试验过程中收集到的齿侧间隙退化数据和定位精度退化数据，考虑到试验成本和机床厂的保密要求，采用块Bootstrap法和插值法对原始数据进行了适当处理，仅保留数据的统计学特性.

由图7和图8可以看出，随着时间的推移，齿侧间隙和定位精度两个性能退化量存在波动和随机性，因此，用随机效应下Wiener过程来表示齿侧间隙的退化量和定位精度的退化量. 齿侧间隙退化有明显的非线性特性，因此认为其性能退化模型中的\Lambda (t)=\Lambda (t,d)=t^d，而定位精度退化轨迹浮动较小且整体退化趋势更符合线性过程. 因此，用随机效应下非线性Wiener过程来描述齿侧间隙退化过程，用随机效应下线性Wiener过程描述定位精度退化过程. 对图7和图8中的样本1、样本2和样本3，设样本序号为i，则i=\mathrm{1,2},\mathrm{3}，Y_{\mathrm{1}}(t_{ij})表示第i个样本在t_j时刻的齿侧间隙退化数据，\mathrm{\Delta }{y}_{{\mathrm{1}}_{ij}}=Y_{\mathrm{1}}(t_{ij})-Y_{\mathrm{1}}(t_{ij-\mathrm{1}})为样本i的齿侧间隙退化增量；Y_{\mathrm{2}}(t_{ij\mathrm{\text{'}}})表示第i个样本在t_{j\mathrm{\text{'}}}时刻的定位精度退化数据，\mathrm{\Delta }{y}_{{\mathrm{2}}_{ij\mathrm{\text{'}}}}=Y_{\mathrm{2}}(t_{ij\mathrm{\text{'}}})-Y_{\mathrm{2}}(t_{ij\mathrm{\text{'}}-\mathrm{1}})为样本i的定位精度退化增量. 为了满足该机床的功能要求，通过分析可得齿侧间隙退化过程的失效阈值为C_{\mathrm{B}}=\mathrm{280}，定位精度退化过程的失效阈值为C_{\mathrm{P}}=\mathrm{500}.

根据观测数据对性能退化模型中的固定参数进行求解. 假定随机效应下非线性Wiener过程的漂移系数\mu是常数，则该随机过程性能退化量的概率密度函数可以表示为：

f(\mathrm{\Delta }{y}_{ij}\left|\mu \right.)={\int }_{\mathrm{0}}^{+\mathrm{\infty }}f(\mathrm{\Delta }{y}_{ij}\left|\mu ,{\sigma }^{\mathrm{2}})\right.g_{\mathrm{1}}({\sigma }^{\mathrm{2}})\mathrm{d}{\sigma }^{\mathrm{2}}=

\frac{\sqrt[]{b}\mathrm{\Gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}-a\right)}{\sqrt[]{\mathrm{2}\mathrm{\pi }\mathrm{\Delta }\Lambda (t)}\mathrm{\Gamma }(a)}{\left\{\left.b^{\mathrm{2}}+\frac{b{\left[\mathrm{\Delta }y-\mu \mathrm{\Delta }\Lambda (t)\right]}^{\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }\Lambda (t)}\right\}\right.}^{a-\mathrm{1}/\mathrm{2}}（33）

固定参数\left\{a\right.,b,d\}以及漂移系数\mu的先验分布为无信息先验分布，采用MCMC-Gibbs抽样算法估计4个未知参数的值，抽取50 000个样本，考虑到前期迭代的不稳定性，设置烧入期为8 000.齿侧间隙和定位精度退化模型固定参数求解迭代过程分别如图9和图10所示.为了减少初始赋值对参数估计结果的影响并更好地判断其收敛性，生成两条MCMC链，给它们赋予不同的初始值，两条不同的MCMC链分别用红线与蓝线表示.

两个性能退化模型中的固定参数估计结果如 表1所示.

对融合失效寿命数据的先验分布超参数进行求解. 元动作单元的可靠性高，失效寿命数据较难获得，用于确定先验分布的蜗杆转动元动作单元失效寿命数据将结合整机故障数据以及相似元动作单元的寿命信息获得. 经分析，蜗杆转动元动作单元的历史失效寿命数据分别为5 150 h、5 469 h、5 038 h. 根据式（11），利用MATLAB中的fminsearch函数求解参数，结果如表2所示.

将\{{\theta }_{{\mu }_{\mathrm{1}}},{\lambda }_{{\mu }_{\mathrm{1}}},{\theta }_{{\mu }_{\mathrm{2}}},{\lambda }_{{\mu }_{\mathrm{2}}}\}的先验分布代入式（20），可以得到融合失效寿命信息的齿侧间隙退化下的后验分布与定位精度退化下的后验分布，利用MCMC-Gibbs抽样法求解这两个后验分布.齿侧间隙和定位精度退化下融合模型未知参数迭代求解过程分别如图11和图12所示，其中红线与蓝线分别表示具有不同初始值的两条MCMC链. 两个退化模型的后验分布参数估计结果如表3所示.

将齿侧间隙退化模型中的未知参数与定位精度退化模型中的未知参数代入式（10），可以得到基于单性能退化的可靠度函数R_{\mathrm{1}}(t)和R_{\mathrm{2}}(t). 根据它们之间的失效相关性进行分析. 图13为失效分布函数F₁和F₂的频率分布直方图.

由前文关于Copula函数的介绍可以得知，椭圆族Copula（包括Gaussian Copula、t-Copula）适用于描述线性相关且上下尾部有对称性的变量关系；阿基米德族Copula（包括Gumbel Copula、Clayton Copula、Frank Copula）能够描述变量间的非对称依赖关系，对于描述尾部相关性较强的数据效果较好. 根据图13的直方图可以排除椭圆族Copula，从阿基米德族Copula中选取合适的Copula函数. 采用极大似然法对Copula函数中的参数进行求解，3种Copula函数参数估计值及AIC、BIC评价准则计算结果如表4所示.

由模型评价准则可以看出，Gumbel Copula的值最小，二维 Gumbel Copula函数的密度函数不具有对称性，上尾部较下尾部高，适用于描述上尾部相关特征的变量，与图13频率分布直方图的观察结果相符，因此本文选用Gumbel Copula函数描述两个失效分布函数之间的相关性.

融合寿命信息的齿侧间隙退化与定位精度退化联合可靠度模型为：

R^C(t)=R_{\mathrm{1}}(t)+R_{\mathrm{2}}(t)-\mathrm{1}+

C\left[F_{T_{\mathrm{1}}}(t),F_{T_{\mathrm{2}}}(t);\theta \right]=R_{\mathrm{1}}(t)+R_{\mathrm{2}}(t)-\mathrm{1}+

\mathrm{e}\mathrm{x}\mathrm{p}\left\{-{\left\{{\left[-\mathrm{l}\mathrm{n}{F}_{T_{\mathrm{1}}}(t)\right]}^{\theta }+{\left[-\mathrm{l}\mathrm{n}{F}_{T_{\mathrm{2}}}(t)\right]}^{\theta }\right\}}^{\mathrm{1}/\theta }\right\}（34）

5.2.2 模型验证与可靠性评估

本文采用交叉验证的方法对所提出的融合寿命信息的性能退化模型进行验证，并利用模型对比的方法验证考虑多性能退化失效相关关系的合理性，在此基础上，使用本文提出的模型完成蜗杆转动元动作单元的可靠性评估. 齿侧间隙和定位精度的前95%退化数据作为模型参数估计过程中的观测数据，后5%的退化数据作为模型验证点对退化模型进行交叉验证. 验证结果分别如图14和图15所示，箱形图的中间线对应预测结果的均值，上下边缘分别对应预测结果90%置信区间上下限.

从图14和图15可以看出，箱形图基本可以涵盖实际观测值，说明预测值与真实值较为接近，本文所提出的性能退化模型和参数估计方法可以很好地描述性能退化过程并且拥有较强的性能退化预测能力.

为了验证本文所提方法的可行性和优越性，采用不同模型求解蜗杆转动元动作单元的可靠度.

建立如下3个不同可靠性模型：

1） 不考虑失效相关性的基于随机效应非线性Wiener过程的可靠性模型M1.

R_{\mathrm{M}\mathrm{1}}(t)=R_{\mathrm{1}}^N(t)\times R_{\mathrm{2}}^N(t)（35）

2） 不考虑失效相关性的基于融合失效寿命数据的随机效应非线性Wiener过程的可靠性模型M2.

R_{\mathrm{M}\mathrm{2}}(t)=R_{\mathrm{1}}(t)\times R_{\mathrm{2}}(t)（36）

式中： R_i （t），i=1，2为两参数威布尔分布的可靠度函数.

3） 考虑失效相关性的基于融合失效寿命数据的随机效应非线性Wiener过程的可靠性模型M3.

R_{\mathrm{M}\mathrm{3}}(t)=R_{\mathrm{1}}(t)+R_{\mathrm{2}}(t)-\mathrm{1}+

\mathrm{e}\mathrm{x}\mathrm{p}\left\{-{\left\{{\left[-\mathrm{l}\mathrm{n}{F}_{T_{\mathrm{1}}}(t)\right]}^{\theta }+{\left[-\mathrm{l}\mathrm{n}{F}_{T_{\mathrm{2}}}(t)\right]}^{\theta }\right\}}^{\mathrm{1}/\theta }\right\}（37）

不同模型下蜗杆转动元动作单元可靠度曲线如图16所示.

由图16可知，模型M1得到的可靠度偏高，主要是因为其性能退化模型与真实情况偏差较大，在此基础上进行可靠性评估准确度会降低. 从试验过程中的性能退化情况来看，M2与M3的可靠度曲线更符合实际情况. 将可靠寿命作为蜗杆转动元动作单元的可靠性指标，当蜗杆转动元动作单元的可靠度达到0.5时所对应的可靠寿命称为中位寿命，蜗杆转动元动作单元的中位寿命为T_{\mathrm{0.5}}=\mathrm{6}\mathrm{792} h. 结合蜗杆转动元动作单元的历史失效寿命进行分析，可以看出，模型M3得到的可靠性评估结果优于模型M2，忽略失效相关关系的可靠性评估方法会低估元动作单元可靠度，影响后续分析. 因此，本文提出的融合寿命信息的性能退化模型更贴合蜗杆转动元动作单元的真实运行状态，考虑多性能退化失效相关性的可靠性模型可以更准确地描述其可靠度变化过程.

6 结 论

本文建立了融合寿命信息与多性能退化信息的可靠性评估模型，主要结论如下：

1） 针对数控机床可靠性建模和分析困难的问题，本文从机床完成其功能的角度出发，依据其运动和动力的传递，将机床分解为一系列元动作单元，准确评估关键元动作单元的可靠性水平对保证机床可靠性具有重要意义.

2） 充分考虑试验样本的个体差异性与共同属性，并且根据实际退化轨迹考虑退化过程存在非线性的情况，使用具有随机效应的非线性Wiener过程建立性能退化模型可以更好地拟合各参数的退化趋势.

3） 针对服从随机效应的非线性Wiener过程的性能退化量，考虑到退化数据不足的情况，提出基于贝叶斯理论的融合失效寿命信息性能退化模型，实例验证表明，该模型可以很好地拟合实际退化过程，能够准确预测未来时间点的退化情况.

4） 考虑多退化失效相关性，根据模型优选准则选用合适的Copula函数进行定量分析，建立的基于Copula函数的多性能退化可靠性模型与模型M1和模型M2进行对比，表明考虑失效相关性的模型具有更高的可靠性评估精度，更符合实际情况. 本文提出的可靠性评估模型准确性更高，能为后续机床可靠性工作提供数据基础和参考依据.

参考文献

[1] 邓超，陶志奎，吴军，等 ．基于性能退化的数控机床剩余寿命预测［J］．机械工程学报，2018，54（17）：181-189．

[2] DENG C， TAO Z K， WU J，et al ．Residual life prediction for NC machine tool based on performance degradation［J］．Journal of Mechanical Engineering，2018，54（17）：181-189．（in Chinese）

[3] 张云，姜楠，王立平 ．基于Wiener过程的数控转台极小子样可靠性分析［J］．清华大学学报（自然科学版），2019，59（2）：91-95．

[4] ZHANG Y， JIANG N， WANG L P. Reliability analysis of NC rotary table based on a Wiener process for extremely small samples［J］. Journal of Tsinghua University （Science and Technology）， 2019， 59（2）： 91-95．（in Chinese）

[5] FENG L， WANG H L， SI X S，et al ．A state-space-based prognostic model for hidden and age-dependent nonlinear degradation process［J］. IEEE Transactions on Automation Science and Engineering， 2013， 10（4）： 1072-1086．

[6] TANG S J， GUO X S， YU C Q，et al ．Accelerated degradation tests modeling based on the nonlinear Wiener process with random effects［J］．Mathematical Problems in Engineering，2014，2014：560726．

[7] 王汉禹， 陈震， 周笛， 等 ．基于核函数-Wiener过程的轧辊非线性退化建模与剩余寿命预测［J］. 上海交通大学学报， 2023，57（8）： 1037-1045．

[8] WANG H Y， CHEN Z， ZHOU D，et al ．Nonlinear degradation modeling and residual life prediction for rollers based on kernel-based Wiener process［J］. Journal of Shanghai Jiao Tong University， 2023， 57（8）： 1037-1045．（in Chinese）

[9] 张忠文， 彭翀， 车众元， 等. 多阶段退化数据融合的伺服驱动单元可靠性建模［J/OL］. 北京航空航天大学学报.（2023-10-11）［2023-12-25］.

[10] ZHANG Z W， PENG C， CHE Z Y， et al. Servo drive unit reliab-ility modeling with multi-stage degradation data fusion［J/OL］. Journal of Beijing University of Aeronautics and Astronautics. （2023-10-11）［2023-12-25］.in Chinese）

[11] 陈传海， 王成功，杨兆军，等 ．数控机床可靠性建模研究现状及发展动态分析［J］．吉林大学学报（工学版），2022，52（2）：253-266．

[12] CHEN C H， WANG C G， YANG Z J，et al ．Research status and development trend analysis of reliability modeling of CNC machine tools［J］．Journal of Jilin University （Engineering and Technology Edition），2022， 52（2）： 253-266．（in Chinese）

[13] 彭卫文 ．重型数控机床可靠性建模与评估技术研究［D］．成都：电子科技大学，2016．

[14] PENG W W ．Research on reliability modeling and evaluation technology of heavy CNC machine tools［D］．Chengdu：University of Electronic Science and Technology of China，2016．（in Chinese）

[15] 杨文. 基于多源数据融合的数控铣床可靠性分析方法研 究［D］．西安： 西安电子科技大学，2021．

[16] YANG W ．Research on reliability analysis method of CNC milling machine based on multi-source data fusion［D］．Xi’an：Xidian University，2021．（in Chinese）

[17] 孙博 ．基于多源信息融合的数控机床关键功能部件可靠性评估［D］．长春：吉林大学，2021．

[18] SUN B ．Reliability evaluation of key functional components of CNC machine tools based on multi-source information fusion［D］．Changchun：Jilin University，2021．（in Chinese）

[19] 李健，慕宗燚，杜彦斌，等 ．基于元动作的数控机床可用性分析方法研究［J］．中国机械工程，2023，34（15）：1848-1855．

[20] LI J， MU Z Y， DU Y B，et al ．Research on availability analysis method of CNC machine tools based on meta-action［J］．China Mechanical Engineering，2023，34（15）：1848-1855．（in Chinese）

[21] 张根保，冉琰，庾辉，等 ．机械产品可靠性研究的新进展：元动作可靠性理论［J］．制造技术与机床，2022（1）：53-59．

[22] ZHANG G B， RAN Y， YU H，et al ．New progress in reliability research of mechanical products—reliability theory of meta- action［J］. Manufacturing Technology & Machine Tool， 2022（1）：53-59．（in Chinese）

[23] YU H， ZHANG G B， RAN Y ．A more reasonable definition of failure mode for mechanical systems using meta-action［J］．IEEE Access， 2019， 7： 4898-4904．

[24] LI X L， RAN Y， ZHANG G B. Optimization of equal-cycle maintenance strategy considering imperfect preventive maintenance［J］. Proceedings of the Institution of Mechanical Engineers，Part E：Journal of Process Mechanical Engineering，2022，236（4）：1392-1402．

[25] 张根保，张定飞，冉琰，等 ．基于Gamma和混合Copula的元动作单元性能可靠性分析［J］．湖南大学学报（自然科学版），2021，48（4）：113-125．

[26] ZHANG G B， ZHANG D F， RAN Y，et al ．Performance reliability analysis of meta-action unit based on Gamma process and hybrid Copula function［J］．Journal of Hunan University （Natural Sciences），2021，48（4）：113-125．（in Chinese）

[27] 鞠萍华，谷豪东，许文林，等 ．考虑外界影响的元动作单元故障概率模型［J］．湖南大学学报（自然科学版），2021，48（10）：21-28．

[28] JU P H， GU H D， XU W L，et al ．Failure probability model of meta-action unit considering external influence［J］．Journal of Hunan University （Natural Sciences），2021，48（10）：21-28．（in Chinese）

[29] ZHU Y， LIU S W， WEI K X，et al ．A novel based-performance degradation Wiener process model for real-time reliability evaluation of lithium-ion battery［J］．Journal of Energy Storage，2022，50：104313．

[30] 董庆来，王伟伟，司书宾 ．基于性能指标融合的随机退化系统竞争失效分析［J］．西北工业大学学报，2021，39（2）：439-447．

[31] DONG Q L， WANG W W， SI S B ．Competitive failure analysis of a stochastic degradation system based on performance characteristics fusion［J］．Journal of Northwestern Polytechnical University， 2021， 39（2）： 439-447．（in Chinese）

[32] LIU M J， DONG Z S， SHI H，et al ．Remaining useful life estimation of fan slewing bearings in nonlinear Wiener process with random covariate effect［J］．Shock and Vibration， 2022，2022：5441760．

[33] TANG S J， XU X D， YU C Q，et al ．Remaining useful life prediction with fusing failure time data and field degradation data with random effects［J］．IEEE Access， 1964， 8： 11964-11978．

[34] 张点，邢云燕，蒋平 ．基于可靠性增长的可靠性鉴定试验方案［J］．系统工程与电子技术， 2023， 45（11）： 3699-3705．

[35] ZHANG D， XING Y Y， JIANG P ．Reliability qualification test plan based on reliability growth［J］．Systems Engineering and Electronics， 2023， 45（11）： 3699-3705．（in Chinese）

[36] 苏春， 汪千程 ．基于改进多层级嵌套Copula方法的机电产品可靠性评估［J］．东南大学学报（自然科学版）， 2022， 52（5）：981-989．

[37] SU C， WANG Q C ．Reliability evaluation for electromechanical products based on improved multi-level nested Copula method［J］.Journal of Southeast University （Natural Science Edition），2022，52（5）：981-989．（in Chinese）

[38] 朱慧婷. 基于Copula相关结构的电力系统负荷预测方法研 究［D］. 北京： 北京交通大学， 2018.

[39] ZHU H T. Research on power system load forecasting based on Copula-related structure［D］. Beijing：Beijing Jiaotong University， 2018. （in Chinese）

[40] ZHANG Z， GAO D， GUAN T， et al. A reliability evaluation method for Gamma processes with multiple random effects［J］. Machines， 2023， 11（9）： 905.

基金资助

国家自然科学基金资助项目（52275473,51835001）