王浩偉，徐廷學，劉 勇

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

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

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

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

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

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

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

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

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

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

2 確定隨機參數(shù)的先驗分布

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

2.1 加速應力下的參數(shù)值估計

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

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

2.2 加速因子估計及參數(shù)值折算

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

周源泉等［9-11］指出，AFij為一個與可靠度無關的常數(shù)，是保證產(chǎn)品在加速試驗中失效機理不變的充分必要條件.對于任意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 為試驗產(chǎn)品總數(shù).

2.3 確定參數(shù)的分布類型

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

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

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

2.4 估計隨機參數(shù)的超參數(shù)值

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

對于LogNorm 分布，似然函數(shù)為

3 Bayesian統(tǒng)計推斷

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

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

可得α＊、β＊的邊緣密度函數(shù)為

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

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 確定隨機參數(shù)的先驗分布

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

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

表1 加速應力下參數(shù)（μ，σ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

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

產(chǎn)品的加速因子表達式如下：

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

通過Anderson-Darling 統(tǒng)計量分別對αh0、βh0備選分布類型進行擬合優(yōu)度檢驗，結(jié)果如表4所示.

表2 AFi0值Tab.2 Values of AFi0

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

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

由表4可知，αh0最優(yōu)服從Gamma分布，βh0最優(yōu)服從LogNorm分布.若使用參數(shù)的共軛先驗分布，則須指定βh0服從Gamma分布.在本實例中，βh0與Gamma分布的擬合效果一般，使用共軛先驗分布進行Bayesian統(tǒng)計推斷不合適，所以采用參數(shù)的非共軛先驗分布.根據(jù)式（17）、（18）可得超參數(shù)的極大似然估計值，參數(shù)的先驗分布確定為α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個測量數(shù)據(jù)（6≤i≤10），預測產(chǎn)品剩余壽命.若僅利用現(xiàn)場測量數(shù)據(jù)，則參數(shù)估計值及剩余壽命預測值如表5所示.若使用Bayesian統(tǒng)計推斷的方法，則可得參數(shù)的后驗均值和剩余壽命預測值，如表6所示，其中通過Bootstrap方法給出95%置信區(qū)間.

表5 僅利用現(xiàn)場測量數(shù)據(jù)所得估計值Tab.5 Estimates only from field measuring data

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

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

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

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

5 結(jié) 論

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

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

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

（

）：

［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］葛蒸蒸，李曉陽，姜同敏.費用約束下含應力優(yōu)化的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］周源泉，翁朝曦，葉喜濤.論加速系數(shù)與失效機理不變的條件（Ⅰ）：壽命型隨機變量的情況［J］.系統(tǒng)工程與電子技術(shù)，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］.推進技術(shù)，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］王浩偉，徐廷學，趙建忠.融合加速退化和現(xiàn)場實測退化數(shù)據(jù)的剩余壽命預測方法［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］王浩偉，徐廷學，周偉.綜合退化數(shù)據(jù)與壽命數(shù)據(jù)的某型電連接器壽命預測方法［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.