基于隨機參數Gamma過程的剩余壽命預測方法

王浩偉，徐廷學，劉 勇

（1.海軍航空工程學院 兵器科學與技術系，山東 煙臺264000；2.海軍航空工程學院 接改裝訓練大隊，山東 煙臺264000）

近些年，越來越多的高可靠性產品出現在軍工、航空、航天等領域，這類產品在投入使用之前通常利用加速試驗預測其可靠性水平及進行定壽.對于某些具有老化特征的高可靠性產品，加速老化試驗是一種效率較高的試驗手段［1-3］，因為這類產品的某些性能指標會隨時間老化，對性能老化量進行測量、建模，則可推測出產品的壽命信息.例如某型導彈電連接器的接觸電阻在使用過程中逐漸增大，當接觸電阻老化量達到失效閾值時發生失效.

在產品投入使用之后，準確預測產品的剩余壽命是實現產品視情維修、提高裝備戰備完好率的關鍵.由于產品個體之間不可避免存在差異，加速老化試驗得到的產品總體信息不能反映出個體差異，無法推測個體的真實剩余壽命.高可靠性產品在額定應力下的老化速率很低，現場測量到的性能老化數據難以對其老化趨勢正確建模，也不能準確地預測個體的剩余壽命.研究如何融合加速老化數據和現場測量數據進行準確的個體剩余壽命預測，對有效地開展視情維修、精確化保障等具有較強的工程指導意義.

產品的老化過程具有不確定性，適合采用隨機過程進行建模，Gamma過程被廣泛用于對具有單調變化趨勢的老化過程進行建模［4-5］.為了描述個體差別，可以使用具有隨機參數的Gamma過程；為了便于Bayesian統計推斷，大多采用參數的共軛先驗分布假定［6-7］：形狀參數為一常量，尺度參數為一服從Gamma分布的隨機變量.這種假定具有很好的統計特性，但難免會出現形狀參數與假定的Gamma分布擬合不理想的情況.本文研究使用參數的非共軛先驗分布，其中尺度參數和形狀參數都設為隨機變量，通過擬合優度檢驗確定它們的最優分布類型.

綜上所述，本文針對老化規律符合Gamma過程的產品，將加速老化數據作為先驗信息，將產品個體在額定應力下的測量數據作為現場信息，使用參數的非共軛先驗分布進行Bayesian統計推斷，進而預測個體剩余壽命.

1 基于Gamma過程的剩余壽命模型

設Y（t）表示產品在時刻t的性能老化量并且有Y（0）＝0，如果Y（t）為Gamma過程，則其獨立、非負增量ΔY（t）＝Y（t＋Δt）-Y（t）服從如下Gamma分布：

式中：Ga（·）表示Gamma分布函數，α（α＞0）為形狀參數，β（β＞0）為尺度參數，Λ（t）為時間t的函數并且有Λ（0）＝0.根據Gamma分布的可加性可知，Y（t）～Ga（α·Λ（t），β），概率密度函數（PDF）為

當Y（t）首次達到失效閾值D 時定義為樣品失效，可靠度函數可由下式得出：

由于式（4）難以進行統計計算，一般采用BS分布對F（t；D）近似擬合［8］.進行如下時間尺度轉換z ＝Λ（t），CDF可以被近似表示為

2 確定隨機參數的先驗分布

由于存在個體差異，每個產品的老化過程并不一致，老化模型的參數值也不同，但對于出自同一總體的個體，可以認為參數值服從同一分布函數.設形狀參數αi和尺度參數βi 分別服從αi～fα（a，b），βi～fβ（c，d），且fα（a，b）與fβ（c，d）之間不相關，π（α，β）為α和β 的聯合概率密度函數，則有π（α，β）＝fα（a，b）·fβ（c，d）.

2.1 加速應力下的參數值估計

設S0為正常工作應力，Sk為第k 個加速應力，yijk為第j 個產品在Sk下的第i 次測量值，tijk為對應的測量時刻，其中i＝1，2，…，n1；j＝1，2，…，n2；k＝1，2，…，n3.Δyijk＝yijk-y（i-1）jk為 老 化 量 的 增量，ΔΛ（tijk）＝Λ（tijk）-Λ（t（i-1）jk）為時間增量.根據每個樣品的n1次測量數據建立似然函數為

將測量數據（Δyijk，ΔΛ（tijk））代入式（8），可以估計每個樣品參數值（αjk，βjk），形式如下：

2.2 加速因子估計及參數值折算

加速因子的一種等效定義可以由Nelson假設推出，設Fi（ti）、Fj（tj）分別為產品在應力Si和Sj下的累積失效概率，如有Fi（ti）＝Fj（tj），則可將Si相當于Sj的加速因子定義為

周源泉等［9-11］指出，AFij為一個與可靠度無關的常數，是保證產品在加速試驗中失效機理不變的充分必要條件.對于任意ti恒有下式成立：

將式（3）代入式（10），并設Λ（t）＝t，可得

若要保證式（11）恒成立，須滿足：

從而可得

可知，β在各應力下不變，而α隨著應力S 的變化而改變.α與S 的關系可以用加速模型來描述，假設應力為溫度T，加速模型為Arrhenius模型，α 可以表示為

將式（14）代入式（13），可得

將加速應力下的αjk代入式（14），通過最小二乘估計可得，進一步可由式（15）計算出Tk相對于T0的加速因子AFk0.由αh0＝αjk／Kk0可得αjk在T0下的折算值αh0，由βh0＝βjk 可得βjk 在T0下的折算值βh0，h＝1，2，…，m，m 為試驗產品總數.

2.3 確定參數的分布類型

使用Anderson-Darling 擬合優度檢驗的方法確定αh0、βh0的最優擬合分布類型.因為Exponential分布、Normal分布、LogNormal分布、Weibull分布、Gamma分布包含了絕大多數參數分布情況，故在這5種分布中選擇最優擬合的分布類型.

Anderson-Darling統計量的AD 值越小，說明備選分布類型與數據擬合得越好［12］，計算式如下：

式中：Fn（x）為經驗分布函數，F（x）為累積分布函數.此外，可以利用Anderson-Darling 對參數的分布類型假設進行檢驗，若Anderson-Darling的p 值大于顯著性水平0.05，則可以接受原假設.

2.4 估計隨機參數的超參數值

假設參數的先驗分布確定為αh0～Ga（a，b），βh0～LN（c，d），其中超參數a、b 為Gamma分布的形狀參數和尺度參數，超參數c、d 為LogNorm 分布的對數均值和對數標準差.根據各自分布的PDF建立極大似然方程.對于Gamma 分布，似然函數［13］為

對于LogNorm 分布，似然函數為

3 Bayesian統計推斷

某個運行中的產品與老化試驗中的樣品出自同一總體，參數值服從同一分布，可以使用Bayesian方法估計參數的后驗均值.

設經過定期檢測得到某一個體的i個性能老化數據為Y＝［Y（t1），Y（t2），…，Y（ti）］，π（α＊，β＊｜Y）為聯合后驗PDF，L（Y｜αh0，βh0）為極大似然函數，則有

可得α＊、β＊的邊緣密度函數為

進一步通過計算α＊、β＊的期望得到隨機參數的后驗均值：

4 實例應用

4.1 加速老化試驗

某型導彈電連接器的主要失效原因是接觸電阻老化，老化速率主要受溫度和濕度的影響，溫濕度可以影響插針表面氧化物的生成，氧化物的堆積導致接觸電阻不斷增大，最終引起接觸失效［15］.試驗中采用恒定應力加速老化試驗方式，將溫度、濕度作為綜合加速應力.

1）試驗樣本量為24，平均分配在3組應力下：S1（t1＝90 ℃，RH1＝75%），S2（t2＝90 ℃，RH2＝95%），S3（t3＝130 ℃，RH3＝95%）.

2）S1下對樣品共進行15次測量，測量間隔為96h；S2下進行12次測量，測量間隔為72h；S3下進行10次測量，測量間隔為48h.

3）在正常工作應力S0（t0＝55℃，RH0＝60%）下對樣品的性能老化量進行測量，該型電連接器接觸電阻的失效閾值為5mΩ.

4.2 確定隨機參數的先驗分布

測量到的樣品老化情況如圖1～3所示，對各樣品是否服從Gamma老化過程進行檢驗.根據文獻［8］可 知，Δyiβ／（Δtiα）近 似 服 從 均 值 為1-1／（9Δtiα）、方 差 為1／（9Δtiα）的 正 態 分 布.將Δyiβ／（Δtiα）正態標準化后，可以通過判斷是否服從標準正態分布得出老化過程是否為Gamma過程的結論.通過Anderson-Darling統計量進行置信度為95%的標準正態分布假設檢驗，證明所有樣品的老化過程都服從Gamma過程.

加速應力下所有樣品參數的極大似然估計值（MLE）如表1所示.

表1 加速應力下參數（μ，σ2）的極大似然估計值Tab.1 MLEs of（μ，σ2）of each sample at accelerated stresses

圖1 S1 下樣品老化曲線Fig.1 Aging curves of samples at S1

圖2 S2 下樣品老化曲線Fig.2 Aging curves of samples at S2

圖3 S3 下樣品老化曲線Fig.3 Aging curves of samples at S3

對于溫度、濕度綜合應力，可以用廣義艾林模型描述α與兩應力之間的關系：

產品的加速因子表達式如下：

利用AFi0對表1中的各參數估計值進行折算，結果如表3所示.

通過Anderson-Darling 統計量分別對αh0、βh0備選分布類型進行擬合優度檢驗，結果如表4所示.

表2 AFi0值Tab.2 Values of AFi0

表3 參數的折算值（αh0，βh0）Tab.3 Transformed（αh0，βh0）

表4 αh0、βh0在各分布下的AD值Tab.4 ADs ofαh0，βh0at different distribution types

由表4可知，αh0最優服從Gamma分布，βh0最優服從LogNorm分布.若使用參數的共軛先驗分布，則須指定βh0服從Gamma分布.在本實例中，βh0與Gamma分布的擬合效果一般，使用共軛先驗分布進行Bayesian統計推斷不合適，所以采用參數的非共軛先驗分布.根據式（17）、（18）可得超參數的極大似然估計值，參數的先驗分布確定為αh0～Ga（15.766，5.887×103），βh0～LN（3.097，0.378）.

4.3 預測剩余壽命

額定應力S0下運行的某電連接器每隔3個月檢測一次，共獲取10組接觸電阻測量值：0.233、0.494、0.731、1.037、1.348、1.512、1.830、2.203、2.529、2.811mΩ.依次使用前i個測量數據（6≤i≤10），預測產品剩余壽命.若僅利用現場測量數據，則參數估計值及剩余壽命預測值如表5所示.若使用Bayesian統計推斷的方法，則可得參數的后驗均值和剩余壽命預測值，如表6所示，其中通過Bootstrap方法給出95%置信區間.

表5 僅利用現場測量數據所得估計值Tab.5 Estimates only from field measuring data

表6 使用Bayesian統計推斷方法所得估計值Tab.6 Estimates from Bayesian method

由表5、6可知，僅利用現場數據時，估計值與預測值的置信區間較大，說明估計精度不高.利用Bayesian統計推斷時，相同估計值的置信區間比僅利用現場數據時小，說明估計精度得到了提高.產品在第i次（6≤i≤10）檢測時，利用Bayesian方法得到剩余壽命的PDF和CDF，如圖4所示.可知，隨著現場數據的增多，PDF曲線和CDF曲線逐漸變窄，說明剩余壽命的預測精度越來越高.圖5給出3種方法的剩余壽命預測曲線：1）僅利用現場數據；2）僅利用先驗數據；3）Bayesian統計推斷.利用方法1）和2）得出的預測曲線之間有較明顯的差異，說明個體和總體的壽命特征不一致，不能相互替代.利用方法3）得出的預測曲線在方法1）和方法2）之間，是因為對現場數據和先驗數據進行了融合，該預測結果的可信度較高.由各預測曲線的間距可知，與先驗數據相比，現場數據對融合預測結果的影響更大.

圖4 剩余壽命概率密度函數變化情況Fig.4 Plot of probability density function of residual life

圖5 利用3種方法得到的剩余壽命預測值Fig.5 Residual life obtained from three different methods

5 結 論

（1）加速老化試驗已成為評估產品壽命的重要途徑，基于Bayesian的剩余壽命預測方法融合了加速老化數據，提高了預測結果的可信度，可以實現預測值的實時更新，為產品的視情維修提供了重要參考.

（2）共軛先驗分布假定只將Gamma過程的尺度參數看作隨機變量，并且要求隨機參數服從Gamma分布，這可能會出現擬合不理想的情況.共軛先驗分布不預先假定隨機參數的分布類型，具有良好的工程應用性.

（3）加速因子是分析、處理加速試驗數據的橋梁和紐帶，可以將加速應力下的試驗信息折算到工作應力下處理，為解決復雜可靠性建模與評估問題提供了有益的思路.

（

）：

［1］NELSON W B.Analysis of performance：degradation data from accelerated tests［J］.IEEE Transactions on Reliability，1981，30（2）：149-155.

［2］林瑞進，陳文華，劉娟，等.航天電連接器加速性能退化試驗可行性研究［J］.工程設計學報，2010，17（4）：318-320.LIN Rui-jin，CHEN Wen-hua，LIU Juan，et al.Research on feasibility of accelerated degradation test for aerospace electrical connector［J］.Journal of Engineering Design，2010，17（4）：318-320.

［3］葛蒸蒸，李曉陽，姜同敏.費用約束下含應力優化的CSADT 設 計［J］.北 京 航 空 航 天 大 學 學 報，2011，37（10）：1277-1278.GE Zheng-zheng，LI Xiao-yang，JIANG Tong-min.Planning of CSADT with stress optimization under cost constraint［J］.Journal of Beijing University of Aeronautics and Astronautics，2011，37（10）：1277-1278.

［4］WANG X.Nonparametric estimation of the shape function in a Gamma process for degradation data［J］.The Canadian Journal of Statistics，2009，37（1）：102-118.

［5］TSAI C C，TSENG S T，BALAKRISHNAN N.Optimal design for degradation tests based on Gamma processes with random effects［J］.IEEE Transactions on Reliability，2012，61（2）：604-613.

［6］WANG X L，JIANG P，GUO B，et al.Real-time reliability evaluation for an individual product based on change-point Gamma and Wiener process［J］.Quality and Reliability Engineering International，2014，30（4）：513-525.

［7］VAN NOORTWIJK J M.A survey of the application of Gamma processes in maintenance［J］.Reliability Engineering and System Safety，2009，94（1）：2-21.

［8］WANG H W，XU T X，MI Q L.Lifetime prediction based on Gamma processes from accelerated degradation data［J］.Chinese Journal of Aeronautics，2015，28（1）：172-179.

［9］周源泉，翁朝曦，葉喜濤.論加速系數與失效機理不變的條件（Ⅰ）：壽命型隨機變量的情況［J］.系統工程與電子技術，1996，18（1）：55-66.ZHOU Yuan-quan，WENG Chao-xi，YE Xi-tao.Study on accelerated factor and condition for constant failure mechanism ［J］.Systems Engineering and Electronics，1996，18（1）：55-66.

［10］楊宇航，周源泉.加速壽命試驗的理論基礎（Ⅰ）［J］.推進技術，2001，22（4）：276-278.YANG Yu-hang， ZHOU Yuan-quan. Theoretical foundation of accelerated life testing（Ⅰ）［J］.Journal of Propulsion Technology，2001，22（4）：276-278.

［11］王浩偉，徐廷學，趙建忠.融合加速退化和現場實測退化數據的剩余壽命預測方法［J］.航空學報，2014，35（12）：3350-3357.WANG Hao-wei，XU Ting-xue，ZHAO Jian-zhong.Residual life prediction method fusing accelerated degradation and field degradation data［J］.Acta Aeronautica et Astronautica Sinica，2014，35（12）：3350-3357.

［12］GRACE A W，WOOD I A.Approximating the tail of the Anderson-Darling distribution［J］.Computational Statistics and Data Analysis，2012，56（12）：4301-4311.

［13］TSENG S T，BALAKRISHNAN N，TSAI C C.Optimal step-stress accelerated degradation test plan for Gamma degradation processes［J］.IEEE Transactions on Reliability，2009，58（4）：611-618.

［14］NTZOUFRAS I.Bayesian modeling using WinBUGS［M］.New York：Wiley，2009：12-55.

［15］王浩偉，徐廷學，周偉.綜合退化數據與壽命數據的某型電連接器壽命預測方法［J］.上海交通大學學報，2014，48（5）：702-706.WANG Hao-wei，XU Ting-xue，ZHOU Wei.Lifetime prediction method for missile electrical connector synthesizing degradation data and lifetime data［J］.Journal of Shanghai Jiao Tong University，2014，48（5）：702-706.