滚动部件是旋转机械中应用最为广泛的机械零件，如轴承、滚轮、曲柄等，其性能与工况的好坏直接影响到与之相联系的转轴以及安装在轴上的齿轮乃至机械设备的性能。滚动部件的缺陷会导致机器剧烈振动和产生噪声，甚至会引起设备的损坏，由于滚动部件是最易损坏的元件之一，因此滚动部件的性能状态监测与寿命预测是保证机械设备正常工作的重要途径之一。

以飞机为例，根据某航空公司统计，在2010至2012年间，B737 NG飞机操纵系统中襟缝翼故障次数达78次，仅次于减速板，远高于其他部件。襟缝翼失效模式大部分为襟、缝翼卡阻，卡阻的原因主要为滚轮滑轨的滚动磨损失效^[1]。以铁路为例，轮轨摩擦与磨损是铁路系统运行过程中最常见的问题^[2]，大秦线重伤轨数每年以20%的速度递增。提高滚动表面的耐磨性能是减少磨损、避免磨损失效的有效手段。由于涂层技术具有在不影响基体材料基本机械性能的条件下，可以提高其耐高温、耐磨和耐蚀等机械物理及化学性能等优点，因此耐高温腐蚀的涂层是保证机械装备正常运行的重要手段，其中WC(碳化钨)涂层以其优异的耐磨性广泛应用在航空航天、机械、石油化工等行业。

目前国内外对于碳化钨涂层的摩擦磨损行为研究较少^[3]，已有的研究主要集中在碳化钨涂层的制备、失效机理以及摩擦磨损性能的影响因素等方面^[4-9]，而在碳化钨涂层的寿命预测方面的研究几乎没有。但国内外在其他涂层的寿命预测方面开展了不少研究，其中周益春等^[10-14]在热障涂层的寿命预测方面，研究和提出了基于连续损伤累计的寿命预测模型、基于特定参量监测的寿命预测模型、基于界面氧化的寿命预测模型等。徐滨士等研究了不同载荷下3Cr13不锈钢涂层磨损寿命及其影响因素等^[15]。

滚动表面在交变载荷的长期作用下会产生磨损、疲劳，从而导致退化失效。在实际的磨损过程中，磨料磨损、粘着磨损、疲劳磨损、腐蚀磨损往往同时存在，相互影响，因此难于采用某一特定磨损公式进行简单描述。目前国内外根据退化机理的不同，针对类似的问题，已经建立了多种退化过程模型。由于Gamma过程是一个独立、增量非负的随机过程，目前已经广泛用于研究单调的退化过程，这一特点符合滚动表面碳化钨涂层的磨损退化特点，为此，采用Gamma过程对滚动表面碳化钨涂层磨损退化失效进行建模，在此基础上，结合碳化钨涂层的实时磨损监测数据，进行碳化钨涂层磨损失效的实时剩余寿命预测研究。

1 基于Gamma过程的碳化钨涂层磨损退化模型 1.1 Gamma过程

碳化钨涂层厚度会随着使用循环的增加而减少，因此将滚动表面碳化钨涂层t时刻的厚度定义为退化量X(t)，X(t)随着使用循环的增加不断增加。通常，当厚度减少到一定值时就认为失效，此时滚动表面碳化钨涂层厚度退化量X(t)=r₀达到最大。因此，滚动表面碳化钨涂层的磨损失效符合Gamma过程的特点。

碳化钨涂层退化量X(t)的概率密度函数为

$Ga\left( xv,u \right)=\frac{{{u}^{v}}}{\Gamma \left( v \right)}{{x}^{v-1}}exp\left( -ux \right){{I}_{(0,\infty )}}\left( x \right)$

(1)

式中:v＞0称为形状参数，u>0称为尺度参数。由于时间t>0，v(t)是一个非递减并且右连续的实值函数，并且v(0) =0。

根据Gamma过程特点可知，对于形状参数为v(t)，尺度参数为u>0的Gamma过程C_f具有如下特性：

(1) X(t=0)=0 的概率为1；

(2) 对于任意的 λ₁=9.47×10^-5 ；

X(τ)-X(t)~Ga(v(τ)-v(t),u)

即X(τ)-X(t)服从Gamma分布Ga(v(τ)-v(t),u)。

(3) 具有独立增量，即对任意的t_n>t_n-1>…>t≥0(n>1) ,x(t_n)-x(t_n-1),…，x(t₂)-x(t₁)相互独立。

由于碳化钨涂层的退化量增量是非负的，所以碳化钨涂层退化是一个平稳的Gamma过程，因此，假定v(t)=ct,c>0。所以，X(t)的概率密度函数可以写成

$\begin{align} & {{f}_{X(t)}}\left( x \right)=Ga\left( x|ct,u \right)= \\ & \frac{{{u}^{ct}}}{\Gamma \left( ct \right)}{{x}^{ct-1}}~exp\left( -ux \right){{I}_{(0,\infty )}}\left( x \right) \\ \end{align}$

(2)

同时可以得到碳化钨涂层退化量X(t)的期望E(X(t))=$\frac{ct}{u}$，方差为Var(X(t))=$\frac{ct}{^{{{u}^{2}}}}$。

1.2 Gamma过程参数估计方法

如果在时刻t_i,i=1,…,n，并且0=t₀＜t₁＜t₂＜…＜t_n，对应的磨损退化数据x_i,i-1,…,n，且0=x₀≤x₁≤x₂≤…≤x_n。考虑磨损的增量δ_i=x_i-x_i-1,i=1,2，…,n的似然函数可以写成

$\begin{align} & L\left( {{\delta }_{1}},\ldots ,{{\delta }_{n}}c,u \right)=\underset{i=1}{\overset{n}{\mathop{\prod }}}\,{{f}_{X\left( {{t}_{i}} \right)-X({{t}_{~}}i-1}})\left( {{\delta }_{i}} \right)= \\ & \underset{i=1}{\overset{n}{\mathop{\prod }}}\,\frac{uc[{{t}^{{}}}_{i}~-{{t}^{{}}}_{i-1}]}{\Gamma (c[{{t}_{i}}-{{t}_{i-1}}])}{{\delta }^{c[{{t}_{i}}-{{t}_{i-1}}]-1}}_{i}exp\left( -u{{\delta }_{i}} \right) \\ \end{align}$

(3)

式(3) 分别对c和u求导，可以得到对于c和u的估计值$\hat{c}$和$\hat{u}$的方程组

$\left\{ \begin{align} & \hat{u}=\frac{\hat{c}{{t}_{n}}}{{{x}_{n}}} \\ & \underset{i=1}{\overset{n}{\mathop{\sum }}}\,[{{t}_{i}}-{{t}_{i-1}}\left] \{\varphi (\hat{c} \right[{{t}_{i}}-{{t}_{i-1}}])-log\delta \}= \\ & {{t}_{n}}log\left( \frac{\hat{c}{{t}_{n}}}{{{x}_{n}}~} \right) \\ \end{align} \right.$

(4)

求解方程组(4) 得到c和u的估计值$\hat{c}$和$\hat{u}$。

也可以根据碳化钨涂层磨损量退化数据x_i,i=1,…,n，先得到各监测时间退化量的均值和方差的估计值$\hat{u}$_i，$\hat{\sigma }$²_i，再根据

$\begin{align} & {{{\hat{u}}}_{i}}=\frac{ct}{u} \\ & {{\sigma }^{2}}_{i}=\frac{ct}{{{u}^{2}}} \\ \end{align}$

(5)

对式(5) 两端求对数后，通过线性回归的方法，解得未知参数c和u的估计值$\hat{c}$和$\hat{u}$。

1.3 基于Gamma过程的碳化钨涂层磨损失效模型

设T_r0表示滚动表面碳化钨涂层厚度退化量X(t)首次到达磨损量阈值r₀的时刻，根据Gamma过程的定义可知滚动表面碳化钨涂层磨损失效时间的分布函数可以写成

$\begin{align} & F\left( t \right)=Pr\left\{ {{T}_{{{r}_{0}}}}\le t \right\}=Pr\left\{ X\left( t \right)\ge {{r}_{0}} \right\}\text{ }= \\ & {{\int }^{\infty }}_{x={{r}_{0}}}{{f}_{X(t)}}\left( x \right)dx= \\ & {{\int }^{\infty }}_{x={{r}_{0}}}\frac{{{u}^{~}}ct}{\Gamma \left( ct \right)}{{x}^{ct-1}}exp\left( -ux \right)dx= \\ & \frac{\Gamma \left( ct,{{r}_{0}}u \right)\text{ }}{\Gamma \left( ct \right)} \\ \end{align}$

(6)

式中，对于给定的a>0,x≥0,有Γ(a,x)=∫_z=x^∞z^a-1exp(-z)dz。由可靠性的定义可以得到滚动表面碳化钨涂层的实时性能可靠性为

$\begin{align} & R\left( t \right)=1-Pr\left\{ {{T}_{{{r}_{0}}}}\le t \right\}=1-Pr\left\{ X\left( t \right)\ge {{r}_{0}} \right\}= \\ & {{\int }^{{{r}_{0}}}}_{0}\frac{{{u}^{ct}}}{\Gamma \left( ct \right)}{{x}^{ct-1}}exp\left( -ux \right)dx= \\ & 1-\frac{\Gamma \left( ct,{{r}_{0}}u \right)}{\Gamma \left( ct \right)~} \\ \end{align}$

(7)

从式(7) 中可以看出，对于不同时刻，如果能确定的参数c和u的估计值和，就可以得到此时的滚动表面碳化钨涂层的可靠性。如果滚动表面碳化钨涂层失效时间的分布函数F(t)可导，可得到滚动表面碳化钨涂层失效时间的概率密度,即

$\begin{align} & f\left( t \right)=\frac{\partial }{\partial t}\left[ \frac{\Gamma \left( v\left( t \right),{{r}_{0}}u \right)}{\Gamma \left( v\left( t \right) \right)} \right]= \\ & \frac{\partial }{\partial \upsilon }{{\left[ \frac{\Gamma \left( \upsilon ,{{r}_{0}}u \right)}{\Gamma (\upsilon )} \right]}_{\upsilon =v(t)}}v\prime \left( t \right) \\ \end{align}$

(8)

如果v(t)可导，则有

$\begin{align} & ~f\left( t \right)=\frac{v\prime \left( t \right)}{\Gamma \left( v\left( t \right) \right)}\cdot \\ & {{\int }^{\infty }}_{{{r}_{0}}u}\{log\left( z \right)-\text{ }\psi \left( v\left( t \right) \right)\}{{z}^{v\left( t \right)-1}}exp\left( -z \right)dz \\ \end{align}$

(9)

式中：$\psi \left( a \right)=\frac{\Gamma \prime \left( a \right)}{\Gamma \left( a \right)}=\frac{\partial log\Gamma \left( a \right)}{\partial a}$,a>0。

一般情况下，很难利用解析方法求得此失效时间的概率密度f(t)，因而可以采用仿真的方法求解。具体步骤为先利用蒙塔卡罗仿真模拟退化量的增长过程来求得滚动表面碳化钨涂层的失效时间，由仿真产生的一组失效数据进一步可求产品失效分布的经验函数，并基于此经验函数进行滚动表面碳化钨涂层失效时间概率密度分布的计算。由于利用Gamma过程对碳化钨涂层磨损失效建模过程中无法得到准确的失效阈值的首达时间分布的解析表达式，为了便于验证方法有效性，本文采用给定可靠度阈值为95%的前提下，反求出此时的剩余寿命。

2 碳化钨涂层磨损剩余寿命预测 2.1 Gamma过程参数的贝叶斯更新方法

在试验初期，由于磨损失效数据量较少，如果按照已有少量磨损数据基础上通过增量似然函数方法或对均值方差做线性回归方法来更新Gamma过程的参数c和u，部件需要求解复杂的非线性方程组，同时由于数据量较少会造成较大的误差。然而可以考虑利用已做完试验的样本数据作为正在进行试验的碳化钨涂层磨损剩余寿命预测的先验经验，通过贝叶斯方法更新Gamma过程的参数c和u，从而提高预测的效率与精度。此外，由于需要在更新过程中形成一个先验链，即现在的后验分布可以作为下一次计算的先验分布，持续进行更新过程，这就需要先验分布和后验分布的形式相同，所以这里定义先验分布与似然函数是共轭的。

根据贝叶斯原理可以得到

$\begin{align} & \pi \left( c,u|{{\delta }_{1}},\ldots ,\delta {{~}_{n}} \right)= \\ & \frac{L({{\delta }_{1}},\ldots ,{{\delta }_{n}}|c,u)\pi \left( c,u \right)}{{{\int }^{\infty }}_{0}{{\int }^{~}}{{\infty }_{0}}L({{\delta }_{1}},\ldots ,{{\delta }_{n}}c,u)\pi \left( c,u \right)dcdu} \\ \end{align}$

(10)

式中:L(δ₁,…,δ_nc,u)表示当给定向量(c,u)时退化增量δ₁,…,δ_n的似然函数;π(c,u)表示向量(c,u)的先验分布;π(c,uδ₁,…,δ_n)表示向量(c,u)的后验分布。根据贝叶斯公式，当磨损监测数据有更新时，可以根据新的退化数据来更新(c,u)的后验分布，从而得到实时的正在试验的滚动表面碳化钨涂层的可靠性评估。对于先验分布π(c,u)，可以考虑先固定c，并假设u的先验分布π(uc)是形状参数为α、尺度参数为β的Gamma分布，由Gamma分布的性质可知：u的后验分布是形状参数为(α+ct_n)、尺度参数为(β+x_n)的Gamma分布。由此，贝叶斯公式(10) 可以改写成

$\begin{align} & ~\pi \left( c,u|{{\delta }_{1}},\ldots ,{{\delta }_{n}} \right)=\pi \left( u|c,{{\delta }_{1}},\ldots ,{{\delta }_{n}} \right)\pi \left( c|{{\delta }_{1}},\ldots ,{{\delta }_{n}} \right) \\ & \propto \underset{i=1}{\overset{n}{\mathop{\prod }}}\,\frac{{{u}^{c[{{t}_{i}}-{{t}_{i-1}}]}}}{\Gamma (c[{{t}_{i}}-{{t}_{~}}_{i-1}])}\text{ }{{\delta }_{i}}^{c[{{t}_{i}}-{{t}_{i-1}}]-1}\cdot \\ & exp(-u{{\delta }_{i}})\times \frac{\beta {{~}^{\alpha }}}{\Gamma \left( \alpha \right)}{{x}^{\alpha -1}}exp\left( -u\beta \right)\pi \left( c \right)~ \\ \end{align}$

(11)

当参数c未知时，u的先验分布π(u|c)可以认为是形状参数为α(c)、尺度参数为β(c)的Gamma分布。由于先验分布和似然函数符合共轭条件，从而u的后验均值可以写成

$E(U|c,{{\delta }_{1}},\ldots {{\delta }_{n}})=\frac{\alpha \left( c \right)+c{{t}_{n}}~}{\beta \left( c \right)+{{x}_{n}}}$

(12)

2.2 碳化钨涂层磨损剩余寿命预测

根据贝叶斯公式更新的c和u的估计值$\hat{c}$和$\hat{u}$，得到滚动表面碳化钨涂层磨损可靠性函数

$R\left( t \right)=1-\frac{\Gamma (\hat{c}t,{{r}_{0}}\hat{u})}{\Gamma \left( \hat{c}t \right)}$

(13)

在得到性能可靠性函数的基础上，对于给定可靠性阈值(如95 %)，可以求出此时滚动表面碳化钨涂层的剩余寿命，从而实现碳化钨涂层磨损失效剩余寿命的预测。

3 实例应用

采用的表面碳化钨涂层WC，弹性模量710 000 N/mm²，泊松比0.21，与之相应的陪试轮滚动体表面的材料是镀硬铬GCr15，弹性模量为208 000 N/mm²，泊松比为0.3，接触应力为1 200 MPa。总共进行了15组试验，分别在2 000，4 000，…，60 000循环时，记录磨损量，磨损量阈值r₀=100 μm，试验中的磨损数据如表 1所示。

由于样本选自不同滚动表面碳化钨的滚轮件满足物理上的独立性，并且初始量满足Gamma过程的定义，然后利用Matlab中的ks检验，对磨损增量数据进行Gamma分布检验。经验证，每一组增量数据都符合Gamma分布，也就意味着滚动表面碳化钨的磨损过程满足Gamma过程的要求。

其中，第一组磨损增量服从Ga(x0.798 47.93) 分布，Gamma对应的拟合结果如图 1所示。

采用同样的方法，可以得到各组磨损数据的参数估计值，如表 2所示。

基于表 2的计算结果，以$\hat{c}$的均值0.795作为c初始值，对于3组检验样本，分别利用贝叶斯方法更新Gamma过程的参数u的估计值$\hat{u}$。表 3给出了8 000循环、16 000循环、30 000循环、40 000循环和50 000循环时的Gamma过程的参数$\hat{u}$的更新结果。

在此基础上，按照式(12) 计算检验样本分别在8 000循环、16 000循环、30 000循环、40 000循环和50 000循环时的可靠度，其中第一组样本的可靠度曲线如图 2所示。

图 2中R₁表示在已知8 000循环磨损监测数据所得到的可靠度曲线，R₂表示在已知16 000循环内磨损监测数据所得到的可靠度曲线，以此类推，R₅是检验样本在50 000循环下的可靠性曲线。从图 2中可以看出，可靠度曲线R₁~R₄在8 000到40 000循环时可靠度较低，这与参考样本一致。同时，随着磨损监测数据的不断累积，尤其当已知50 000循环的磨损监测数据时，可靠度曲线已靠近检验样本的真实值，从而实现了对于碳化钨涂层可靠性的实时评估。

根据图 2得到的实时可靠度曲线，在给定可靠度阈值为95%的条件下，就可以得到碳化钨涂层不同磨损循环下的剩余寿命预测值，如表 4所示。为了说明Gamma过程参数的贝叶斯更新方法的有效性，表 4中同时给出了3组检验样本Gamma过程参数未更新下(固定参数)的碳化钨涂层剩余寿命预测值以及实际的剩余寿命值。

从表 4中可以看出，随着磨损监测数据的积累，参数更新下的预测方法预测精度在逐渐提高，相比固定参数的相对误差较小。但是必须指出的是，本文采用给定可靠度阈值为95%的前提下，反求出此时的剩余寿命，从而此时的剩余寿命是在一定安全条件下的估计，与实际存在着一定的误差。

4 结束语

根据滚动表面碳化钨涂层的磨损失效特点，建立了基于Gamma过程的碳化钨涂层磨损失效模型，应用贝叶斯更新方法融合新增的磨损监测数据，对Gamma过程参数进行更新，从而实现了对碳化钨涂层磨损失效的剩余寿命预测。在预测精度上，本文所提的预测方法比之前参数不更新的预测方法有了较大提高，并且剩余寿命预测的精度随着磨损监测数据的增加在不断提高，可以较好地满足工程实际需要。