摘要：动量轮是卫星姿态控制系统的关键部件，对卫星的可靠性及整体寿命有直接影响。针对卫星动量轮的剩余寿命预测问题，提出一种基于Copula函数的多退化量下的寿命预测方法。首先，分析了影响动量轮寿命的关键因素，选择润滑剂剩余量和电流作为退化量；其次，分别对单个退化量进行退化建模，得到动量轮剩余寿命的边缘分布函数；然后，通过Copula函数族来描述多退化量之间的相关性，并对边缘分布进行融合，得到动量轮剩余寿命的联合分布函数；最后，提出基于赤池信息准则（Akaike information criterion，AIC）模型评价的Copula函数选择方法。
关键词：寿命预测；可靠性；性能退化；Copula函数；动量轮
中图分类号：TB114.3 文献标志码：A 文章编号：1005-2615（2015）03-0360-07
Remaining Lifetime Prediction for Momentum Wheel Based on Multiple Degradation Parameters
Liu Shengnan， Lu Ningyun，Cheng Yuehua，Jiang Bin，Xing Yan
Abstract：Momentum wheel （MW） is a critical component in satellite attitude control systems. Its remaining lifetime has a direct impact on the reliability and the useful life of an on-orbit satellite. A lifetime prediction method based on multiple degradation parameters is proposed. Key factors connecting to the MW's lifetime are analyzed. The lubricant remaining amount and current are selected as two degradation parameters of an in-service MW. Then， two marginal distribution functions of MW's remaining lifetime are obtained using the two degradation parameters respectively. After that， Copula function is adopted to describe the correlation between the two degradation parameters. The joint distribution function is then derived by fusing the two marginal distribution functions. Thus， a remaining lifetime prediction model is obtained. Finally， the Akaike information criterion（AIC） is used to select the best Copula function.
Key words： lifetime prediction； reliability； degradation； Copula function； momentum wheel
动量轮是卫星三轴稳定控制系统中的主要执行部件，一般包括轴承组件、电机组件、壳体组件和轮体组件4个部分。动量轮的寿命直接影响卫星的使用寿命，实现动量轮剩余寿命的准确预测，对提高卫星可靠性进而延长卫星使用寿命具有重要意义。
基金项目：国家自然科学基金（61203091，61273171，61374141）资助项目。
收稿日期：2015-02-01；修订日期：2015-04-28
由于动量轮失效机理复杂、价格昂贵，难以进行大样本加速寿命试验，因此目前并没有成熟的动量轮失效机理模型。针对这一类对象，利用没备退化数据建立其性能退化模型是目前被广泛使用的寿命预测方法[1]。现有研究表明，卫星姿态控制系统中，动量轮失效与多种因素相关[2]。外部环境季节性变化、轴承温度、润滑剂系统失效、微振动、电流、转速等冈素都可能影响动量轮寿命[3]。现有的基于数据的剩余寿命预测方法，大多仅考虑单一退化量对剩余寿命的影响，而实际上，动量轮寿命受多个因素的共同影响。譬如，文献[4]研究了润滑剂失效物理特性，建立了润滑剂耗损模型用来预测JB-3动量轮轴承组件寿命；文献[5]以动量轮轴承温度作为动量轮的寿命特征量，利用经验模式分解提取时间变化趋势，建立了动量轮的退化模型，进行在轨动量轮寿命预测；文献[6]使用基于RVMPF的预测方法，以动量轮轴承温度等指标作为动量轮的寿命特征值，建立动量轮的退化模型并进行其寿命预测；文献[7]分析了某型号5个动量轮轴承温度的遥测数据，对其轴温进行性能退化趋势建模，用来预测动量轮剩余寿命；文献[8]指出动量轮失效的关键是供油系统的剩余油量，对动量轮润滑的失效机理进行研究，得到了卫星动量轮的寿命可靠度曲线图；文献[9]研究了动量轮轴系试验数据，发现电流数据的缓慢上升趋势系由轴承摩擦力增大所敛，建立了动量轮电流的时间序列模型，根据其性能退化特性，推导出动量轮寿命预测曲线；文献[10]针对两个相似型号的动量轮产品数据，假设其寿命服从Weibull分布，估计Weibull分布的形状参数，描述了不同型号相似产品之间寿命分布的相似程度；文献[11]给出了在高温边缘润滑状态、季节性变化、温度、微振动等条件下动量轮寿命测试的试验结果；文献[12]建立了轴温、电流和动量轮零部件故障模式之间的网络拓扑结构，利用贝叶斯网络对动量轮进行可靠性建模。 由于目前动量轮失效机理模型尚不明确，很难从小提取某一个退化量作为表征动量轮寿命的决定性因素。因此，本文提出一种基于多退化量的动量轮剩余寿命预测方法。首先，分析了影响动量轮寿命的关键因素，选择润滑剂剩余量和电流作为退化量；其次，对各单退化变量分别建模，得到动量轮剩余寿命的边缘分布函数；然后，通过Copula函数族来描述多退化量问的相关性，得到动量轮剩余寿命的联合分布函数；最后，通过赤池信息准则（Akaike information criterion，AIC）模型评价方法得到最适合的Copula函数；确定最终的多退化量剩余寿命预测模型。仿真实验结果表明了本文方法的可行性和合理性。
1 影响动量轮寿命的关键因素分析
工程实践及大量试验表明，在轨动量轮失效主要原因是其润滑系统失效和轴承失效。润滑系统失效主要有3个影响因素：金属磨损、保持架磨损、润滑剂缺失。对于轴承，轴承几何结构、材料选择、环境因素等都会成为轴承失效的原因[13]。文献[14，15]中给出了轴承失效模式及动量轮主要故障模式。轴承过度供油和过度贫油的情况下.轴温会持续偏高、电流遥测值增大、转速下降；电机驱动电流相短路或断路的情况下，电流遥测值有波动，转速下降，控制指令与电流值关系紊乱。经分析，轴温、电流、转速等性能数据与动量轮故障模式有着复杂的对应关系。文献[11，l6]给出了转速、电流以及温度的监测结果，并发现二者之间存在密切关系，当其中一变量发生变化时，其余两个量也随之发生变化。
考虑到：（1）现有文献已对润滑剂剩余量与动量轮剩余寿命之间的关系进行了深入研究；（2）润滑剂剩余量与轴温、电流与转速之问有清晰的对应关系，但如何描述它们之间的关系尚无研究成果。本文选择润滑剂剩余量和电流作为退化变量，建立基于多退化量的动量轮剩余寿命预测模型。
2 基于多退化量的动量轮剩余寿命预测方法
2.1 单退化量下的寿命预测模型
对于润滑剂，采用文献[5]中提供的润滑剂剩余量的数据；对于电流，根据现有文献[17]，可假定其服从指数分布，仿真生成其性能退化的测试数据。基于所获得的数据，分别建立模型来描述单退化量的退化轨迹，得到动量轮寿命的边缘概率分布，具体步骤如下。
（1）定义动量轮的性能退化函数yi（t）（其中，i=1表示润滑剂剩余量；i=2表示电流）。性能退化函数是指将不同样本在同一时刻所测得的数据看作服从某一分布的随机变量，从退化量分布的角度来描述产品的性能退化规律。由于润滑剂剩余量及电流呈现近似线性衰退，故采用线性退化函数来描述退化过程
式中：参数α、β分别代表退化量的初始值及耗损率，它们相互独立，服从高斯分布，即。显然，yi（t）也服从高斯分布。
（2）由于参数α、β相互独立，所以退化量yi（t）的均值和均方差可以表示为
式（2）中的参数μαi，μβi，σαi，σβi可根据多个动量轮样本的退化量在不同时刻的测量值采用线性回归方法估计得到，详见第3节。
（3）定义失效阈值Di，计算动量轮的寿命失效分布函数Fi（t）。润滑剂剩余量的失效阈值D1可借鉴文献[5]中的取值。假设D1服从正态分布，取值位于区间[Da，Db]，Da、Db为区间边界。电流的失效阈值D2为固定数值，可根据专家经验给出。
基于润滑剂剩余量的动量轮寿命失效分布函数F1 （t）为
式中：P（y1（t）≤D1）表示润滑剂剩余量y1 （t）小于失效阈值D1的概率；Φ表示正态分布函数。
式中：x为数学意义上的广义积分项；因此，式（4）可以表示为
同理，基于电流的动量轮寿命失效分布函数F2（t）为
（4）求动量轮剩余寿命的边缘分布函数Ri（t）。当产品开始使用时，产品的剩余寿命概率为1，即Ri（t）=1（在零时刻，产品总能正常工作）。随着使用时间的增加，产品剩余寿命的概率越来越低，而产品的失效分布函数Fi（t）会随时间而逐渐上升[16]，满足关系：Ri（t）+Fi（t）=1。
因此，基于润滑剂剩余量的动量轮寿命的边缘分布函数为同理，基于电流的动量轮寿命的边缘分布函数为
2.2 基于Copula函数的动量轮剩余寿命预测模型
2.2.1 Copula函数
若设备的各个单退化量相互独立，根据概率论中相互独立的随机变量理论，令单退化量下得到的设备剩余寿命的边缘分布函数相乘，可得到设备剩余寿命的联合分布函数。若设备的各个单退化量相关性已知，即相关系数及协方差已知，可以通过数理统计中的协方差矩阵得到设备剩余寿命的联合分布。若设备的各个单退化量相关性未知，则以上两种做法都不可行[15]。在工程实践中，退化变量之间的相关性往往未知，在给定各单退化量的边缘分布情况下，如何确定它们的联合分布便成了一个非常重要的问题。
卫星动量轮中的润滑剂剩余量与电流之间存在相关性，但相关系数未知。针对这种情况，本文选用Copula函数推导设备剩余寿命的联合分布函数，并提出利用AIC准则，选择最优的Copula函数。表1给出了3种最常用的二元Copula函数，其中，u，v分别是润滑剂剩余量和电流对应的动量轮剩余寿命边缘分布函数，即u=F1（t），v=F2（t）；γ为相关系数，γ取下边界值时，随机变量u，v趋向独立，γ取上边界值时，随机变量u，v趋向于完全相关。
上述3个Copula函数中，不同参数γ下，尾部的非对称性不明显时，或者有多个Copula函数都符合这种尾部特性时，便需要利用模型评价的方法来选择最优的Copula函数模型。
2.2.2 基于Copula函数的寿命预测算法步骤
确定了各单退化量下动量轮剩余寿命的边缘分布后，利用Copula函数计算动量轮剩余寿命的联合分布函数，具体算法步骤如下。 （1）已知润滑剂剩余量的寿命失效边缘分布F1（t）和电流的寿命失效边缘分布F2 （t），令U=F1（t），V=F2（t）。已知U，V取值范围为[0，1]，对它们进行离散化得到（um，vm）（m=1，2，3…，n），绘制其二元直方图，根据二元直方图的形状初步筛选合适的Copula函数。
（2）Copula函数中含有未知参数γ，因此需要进行参数估计。由于润滑剂剩余量及电流的边缘分布已知，本文采用规范化的极大似然估计法（Canonical maximum likehood method，CML）估计参数γ
（3）分别计算Copula密度函数cγ（u，v）和Copula分布函数Cγ（u，v）。
（4）根据概率论中的联合概率分布计算公式，计算基于Copula分布函数的动量轮剩余寿命概率函数R（t）
式中：F1（t）和F2（t）分别为基于润滑剂剩余量和基于电流的动量轮失效分布函数；R1（t）和R2（t）为相应的动量轮剩余寿命概率分布函数；Cγ（u，v）为Copula联合分布函数。
2.2.3 基于AIC的Copula函数选择
实际应用中，由于Copula函数种类较多，对于同一个问题选用不同的Copula函数，可能会导致不同的分析结果。目前还没有通用的方法能选择最合适的Copula函数，来建立联合分布函数，即最符合边缘分布相关性的Copula函数。只能经过初步筛选，选择符合尾部相关性的几个Copula函数。因此，本文提出利用模型评价的方法来确定最优的Copula函数。AIC准则是最常用的模型评价准则，具有良好的适用性和实用性，因此本文选用AIC准则来确定最优的Copula函数。
AIC准则如下
式中：K为被估计参数个数；n为观测样本数目；L为Copula函数的似然值。一般而言，K越小意味着模型越简洁，L越大意味着模型越精确。这里AIC反映了Copula函数拟合原始数据的情况，AIC越小，说明所使用的模型能更好拟合原始数据的联合分布特性。
本文所考虑的3种Copula函数的似然函数如下所示：
Gumbel函数的对数似然函数Frank函数的对数似然函数Clayton函数的对数似然函数
3 仿真验证及分析
本文采用文献[5]中提供的润滑剂剩余量的数据，建立其退化轨迹模型。该数据由地面试验模拟卫星空中运行状态，并在一段时间之后通过解剖动量轮、红外线光谱测量等方法得到。表2给出了轴承润滑剂系统的质量变化数据，包括润滑剂的质量初始值、第8个月的测量值及第11个月的测量值。
对每一时刻的轴承润滑剂系统的均值及方差进行线性拟合，则可以估计出式（2～6）中所用到的相关参数。参数结果为
μα1一5.232，σαl=0.024 38，σβ1=0.412 9，
μβ1=-0.042 19
D1所在的区间为[0.052 3，0.078 5]，得到润滑剂剩余量的寿命分布概率曲线如图1所示。
由于指数分布适用于固有恒定失效率的部件、经老炼试验并进行定期维修的部件[16]，而动量轮电机是具有此特征的部件，所以本文利用指数分布函数仿真生成电流退化数据，如表3所示。
图2给出了基于电流的动量轮剩余寿命概率分布曲线。
将润滑剂剩余量及电流边缘分布的数值范围等分为30个区间，计算边缘分布落在该区间内的频数，得到润滑剂剩余量及电流的频数直方图如图3所示。观察图3可知，该二元频数直方图具有非对称尾部，联合分布密度函数（即Copula函数）应该也具有非对称尾部，因此可以选取二元Gumbel Copula函数或Clayton Copula函数来描述原始数据的相关结构。图中下尾高、上尾低，较符合Clayton Copula函数的特征，但上尾部又并非完全不相关，所以不能排除其符合Gumbel Copula函数和Frank Copula函数的特征。
采用这3种Copula函数所得到的动量轮剩余寿命预测函数的联合分布如图4所示。
从图中可以看出，Gumbel函数和Clayton函数下的动量轮寿命联合分布曲线表现出的曲线特性相似。两者在使用初期（60个月之内）曲线值接近基于润滑剂剩余量的动量轮寿命曲线，之后剩余寿命的概率迅速下降，在120个月之后与基于电流的动量轮寿命曲线基本重合。两者之间不同之处在于在60个月到120个月内Gumbel函数下的动量轮寿命预测概率稍微高于Clayton函数，即Gumbel函数下的动量轮寿命曲线和Clayton相比更接近基于电流的动量轮寿命曲线。Frank函数下的动量轮寿命曲线在60个月寿命概率稍高于基于润滑剂剩余量的动量轮寿命，之后概率迅速降低，在120个月左右寿命概率低于0，这是不符合实际工程应用的，因此得出结论Frank函数不适用于本文所研究的情况，同时也验证了利用二元频数直方图所筛选的适用Copula函数（Gumbel函数和Clayton函数）是可行的。
由图4的结果难以判断哪一种Copula函数所得到的结果更好，因此，必须通过AIC准则对这几种模型进行检验。先用最大似然函数估计法分别得到3种Copula函数的参数值γ，然后计算其似然值，再根据式（12）得到AIC值，判断最优拟合模型，结果见表4。
通过表4比较可以看出在各时刻Clayton函数的AIC最小，说明Clayton Copula函数是用来描述润滑剂剩余量及电流这两个退化变量之间相关性的较好选择，即认为是联合分布里的最优模型。
比较图4（b）中的Clayton Copula联合分布曲线与两个边缘分布曲线发现，对于同一研究对象动量轮来讲，基于润滑剂的动量轮寿命预测曲线符合机械部件的寿命退化情况，而基于电流的动量轮寿命预测曲线符合电子器件的寿命退化情况。但足，两者预测的动量轮寿命有一定差异，因此，有必要对动量轮建立多退化量的联合分布函数。图中基于Clayton Copula函数的动量轮寿命预测曲线，在60个月之前与基于润滑剂的动量轮寿命预测曲线相吻合，而之后与基于电流的动量轮寿命预测曲线相接近。这种情况说明，在动量轮运行初期，寿命退化与润滑剂有很大关系，随着时间变化，动量轮寿命退化与电子元器件的寿命退化更相关。 4 结束语
本文采用数据驱动的方法来研究动量轮寿命预测问题，分析了影响动量轮寿命的主要因素，建立了单退化量及多退化量的动量轮寿命预测模型。本文研究方法与现有的方法相比，主要创新性在于提出一种基于多退化变量的动量轮剩余寿命模型，在两个退化量相关性未知的情况下，采用Copula函数融合两个边缘分布得到动量轮剩余寿命的联合分布，并通过AIC模型评价，选择最合适的Copula函数。结果表明，通过线性退化建模方法得到润滑剂、转速的边缘分布，再利用Copula函数计算联合分布，最终得到的动量轮剩余寿命模型符合实际工程经验。这对于存在多个退化量，且变量之间相关性未知的设备剩余寿命预测具有一定的参考价值。
由于缺乏同一型号同样使用条件下的动量轮寿命数据，尚无法完成预测方法的实证研究。单退化量建模过程中，利用退化量数据的近似线性关系，采用了线性的退化轨迹函数。但在实际应用中，退化轨迹存在非线性，因此本文后续工作还需要研究非线性的退化轨迹函数，使其更符合动量轮退化数据的实际特性。
参考文献：
[1] Ramasso E，Gouriveau R. Remaining useful life estimation by classification of predictions based on a neuro-fuzzy system and theory of belief functions[J].IEEE Trans on Reliability， 2014， 63（2）：555-566.
[2] 祁海铭，程月华，姜斌，等，基于多状态故障的卫星姿态控制系统剩余寿命预测方法[J].南京航空航天大学学报，2015， 47（1）：29-36.
Qi Haiming， Cheng Yuehua， Jiang Bin， et al. Remaining lifetime prediction based on multiple fault states for satellite attitude control system[J]. Journal of Nanjing University of Aeronautics & Astronautics， 2015， 47（1）：29-36.
[3] 康国华，夏青，成婧.基于SoPC的微小卫星姿轨控计算机设计与实现[J].南京航空航天大学学报，2013， 45（6）：763-768.
Kang Guohua，Xia Qing，Cheng Jing. Design of attitude and orbit control computer of micro satellite based on SoPC[J]. Journal of Nanjing University of Aeronautics & Astronautics， 2013，45（6）： 763-768.
[4] Jin G， Liu Q， Zhou J， et al. RePofe： Reliability physics of failure estimation based on stochastic performance degradation for the momentum wheel[J].Engineering Failure Analysis， 2012， 22（2）： 50-63.
[5] Jin G，Matthews D， Fan Y， et al. Physics of failure based degradation modeling and lifetime prediction of the momentum wheel in a dynamic covariate environment[J]. Engineering Failure Analysis， 2013， 28（3）：222-240.
[6] Pecht M， Gu J. Physics-of-failure-based prognostics for electronic products[J]. Transactions of the Institute of Measurement and Control， 2009， 31（3/4）：309-322.
[7] 厉海涛，金光.动量轮维纳过程退化建模与寿命预测[J].航空动力学报，2011，26（3）：622-628.
Li Haitao， Jin Guang. Momentum wheel wiener process degradation modeling and life prediction[J].Journal of Aerospace Power，2011，26（3）：622-628.
[8] 刘强，周经伦，金光，等.基于随机阈值的Gauss＼Brown失效物理模型的动量轮可靠性评估[J].宇航学报，2009， 30（5）：2109-2115.
Liu Qiang， Zhou Jinglun， Jin Guang，et al. Based on random threshold Gauss＼rown failure physical model of reliability assessment of momentum wheel [J].Journal of Astronautics， 2009， 30（5）： 2109-2115.
[9] 刘强.卫星动量轮性能可靠性建模与评估方法研究[D].长沙：国防科学技术大学，2006：14-29.
Liu Qiang. Reliability modeling and evaluation method of satellite momentum wheel's performance[D].Changsha：National University of Defense Technology，2006：14-29. [10] 金光，冯静.长寿命卫星活动部件Bayes-Weibull可靠性评估方法[J].系统工程与电子技术，2009， 31 （8）：2020-2083.
Jin Guang， Feng Jing. Bayes-Weibull reliability assessment method for long life satellite moving component [J]. System Engineering and Electronics，2009，31（8）：2020-2083.
[11] McMahon P， Lavèn R. Results from 10 years of reaction/momentum wheel life testing[C]//11th European Space Mechanisms and Tribology Symposium. [S.l.]：ESMATS，2005： 299-305.
[12] 厉海涛.基于贝叶斯网络的动量轮可靠建模与分析[D].长沙：国防科技大学，2007：35-37.
Li Haitao. Reliability modeling and analysis of momentum wheel based on baysian network[D]. Changsha： National University of Defense Technology，2007：35-37.
[13] 刘良勇.动量轮轴承磨损寿命研究[D].洛阳：河南科技大学，2011：35-36.
Liu Liangyong. Study on wear life of momentum wheel bearing[D]. Luoyang： Henan University of Science and Technology， 2011： 35-36.
[14] 刘文静，刘成瑞.定量与定性相结合的动量轮故障可诊断性评价[J].中国空间科学技术，2011，8（4）：55-63.
Liu Wenjing， Liu Chengrui. The diagnosed evaluation of momentum wheel failure combination of quantitative and qualitative[J]. Chinese Space Science and Technology，2011，8（4）：55-63.
[15] Pecht M， Gu J. Physics-of-failure-based prognostics for electronic products[J]. Transactions of the Institute of Measurement and Control， 2009， 31（3/4）：309-322.[16] 张建勋，胡昌华.多退化变量下基于Copula函数的陀螺仪剩余寿命预测方法[J].航空学报，2014，35（4）：1111-1121.
[16]Zhang Jianxun， Hu Changhua. Multiple degradation variables modeling for remaining usefull life estimation of gyro based on Copula function[J]. Acta Aeronautica et Astronautica Sinica， 2014， 35（4）：1111-1121.
[17] 赵宇.可靠性数据分析[M].北京：国防工业出版社，2011： 21-34.
Zhao Yu. Data analysis and reliability[M]. Beijing： Defense Industry Press，2011： 21-34.