曹宇宇，唐家銀，何平，曹玲

(西南交通大學(xué)數(shù)學(xué)學(xué)院, 四川 成都 611756)

0 引言

隨著科技化水平的不斷進步，很多領(lǐng)域都出現(xiàn)了高可靠性、長壽命的產(chǎn)品，如高性能的電子產(chǎn)品、航天產(chǎn)品，因此用傳統(tǒng)的可靠性試驗會花費大量的時間和成本，且有些傳統(tǒng)的可靠性試驗方法在工程中較難實現(xiàn).為了縮短試驗周期、降低試驗成本，新的可靠性評估方法應(yīng)運而生.其中，加速退化試驗(accelerated degradation testing, ADT)和加速壽命試驗(accelerated life testing, ALT)目前已成為可靠性試驗領(lǐng)域的兩大重要組成部分[1]，受到國內(nèi)外研究學(xué)者的關(guān)注.

ADT是指通過提高應(yīng)力水平來加速產(chǎn)品退化，通過高應(yīng)力下的數(shù)據(jù)來估計產(chǎn)品的可靠性并預(yù)測常應(yīng)力下產(chǎn)品的壽命[2].傳統(tǒng)的加速退化試驗統(tǒng)計分析中通常假設(shè)產(chǎn)品僅有單一失效模式，但在實際工程中，產(chǎn)品存在多種失效模式的可能性很大，任何一個失效模式發(fā)生，產(chǎn)品則發(fā)生失效.從失效模式來看，產(chǎn)品的失效可分為退化失效與突發(fā)失效.在實際工作中，產(chǎn)品發(fā)生失效是由最早出現(xiàn)的失效模式導(dǎo)致的.因此在競爭失效場合下，傳統(tǒng)單一失效模式的統(tǒng)計分析方法并不適用.

關(guān)于競爭失效的研究，最早是由Nelson[3]建立競爭失效模型，通過分析給出了對數(shù)正態(tài)分布下的極大似然估計(MLE).趙建印[4]通過對金屬化膜脈沖電容器的研究，得到了競爭失效下的產(chǎn)品可靠性模型.王華偉等[5]利用貝葉斯模型建立航空發(fā)動機性能退化模型,實現(xiàn)了在不同退化情況下航空發(fā)動機剩余壽命的預(yù)測.安宗文等[6]通過Wiener過程對連續(xù)退化進行刻畫,推導(dǎo)出產(chǎn)品在失效相關(guān)模式下的競爭失效可靠性分析模型.龍哲等[7]通過對退化失效采用Wiener過程來描述產(chǎn)品的性能退化軌跡.上述基于偽失效壽命研究中使用的退化模型大多假設(shè)為服從簡單的一般線性過程或者Wiener過程，并不能很好地擬合退化軌跡，而且基于退化量分布更加關(guān)注的是失效過程本身，比起偽壽命結(jié)果可信度更高.蘇春等[8]則基于性能退化數(shù)據(jù)對競爭失效模型進行了研究和分析，但研究過于簡單，不能夠全面地反映競爭失效問題.另外，對競爭失效下恒定應(yīng)力加速退化試驗也鮮少有較為詳盡的研究.

為了對突發(fā)失效模型進行更好地刻畫，引入基于比例危險的突發(fā)失效模型，比例危險模型最早由COX[9]提出并應(yīng)用于壽命數(shù)據(jù)分析，該模型主要采用非參數(shù)統(tǒng)計方法獲取產(chǎn)品的可靠度估計，EGHBALI[10]改進了COX的比例危險模型，得到了一種融合加速退化因子的比例危險退化模型；蔡忠義等[11]將監(jiān)測時刻作為協(xié)變量，建立了比例危險退化模型.王智明等[12]采用混合Weibull分布并結(jié)合比例危險模型對刀具故障數(shù)據(jù)進行建模.以上方法均是基于退化數(shù)據(jù)，而本研究則對比例危險模型進行了改進，使其應(yīng)用于突發(fā)失效數(shù)據(jù)，并由此獲得基于比例危險的突發(fā)失效模型.

筆者針對競爭失效最一般的形式(產(chǎn)品具有單一退化失效模式和單一突發(fā)失效模式)，對其恒定應(yīng)力加速退化試驗進行建模和分析.首先給出競爭失效場合恒定應(yīng)力加速試驗的突發(fā)與退化競爭模型，在此基礎(chǔ)上分別給出突發(fā)失效模型和加速退化模型，通過對參數(shù)進行估計，實現(xiàn)對競爭失效下恒定應(yīng)力加速退化試驗的統(tǒng)計分析.最后給出應(yīng)用算例.

1 競爭失效場合恒定應(yīng)力的可靠性評估模型

1.1 模型假設(shè)與符號競爭失效場合恒加試驗統(tǒng)計分析基于以下假定：

假設(shè)1：產(chǎn)品具有單一突發(fā)失效模式和單一退化失效模式，其相應(yīng)的突發(fā)失效時間為Tr，退化失效時間為Td，兩者相互競爭；

假設(shè)2：不同時刻、不同應(yīng)力下，產(chǎn)品性能退化量X(t)服從相同的分布族，參數(shù)為時間t和應(yīng)力水平S的函數(shù)；

假設(shè)3：失效閾值D是恒定的；

假設(shè)4：性能退化程度可能會影響突發(fā)失效，也可能不會.

1.2 突發(fā)與退化競爭模型設(shè)在應(yīng)力水平Si(i=1,2,…,k)下，時刻t產(chǎn)品的退化量為X(t)，且當(dāng)X(t)首次達到失效閾值D時產(chǎn)品發(fā)生退化失效.進一步假設(shè)產(chǎn)品突發(fā)失效后功能立刻完全喪失，因此不能繼續(xù)測量其退化量.記退化失效時間為Td,突發(fā)失效時間為Tr，則產(chǎn)品失效時間T為

T=min{Td,Tr}

(1)

時間t內(nèi)產(chǎn)品不失效的概率，即可靠度為

R(t)=P{T>t}=P{Td>t,Tr>t}

(2)

在退化失效模式下，產(chǎn)品的可靠度函數(shù)為

(3)

式中：Rd(t,S)是僅考慮應(yīng)力S下產(chǎn)品退化失效模式時的可靠度，gd(t,xt,S)為應(yīng)力S下退化量在t時刻的分布密度函數(shù)，xt表示t時刻X(t)的取值，即X(t)=xt.

為了得到恒定應(yīng)力下系統(tǒng)競爭失效的可靠度，需要考慮性能退化和突發(fā)失效的相關(guān)性.若突發(fā)與退化相關(guān)，則產(chǎn)品突發(fā)失效的失效率λr(t,xt,S)是時間t、退化量xt和應(yīng)力水平S的函數(shù)，Tr為突發(fā)失效的時間，則突發(fā)失效模式下的可靠度函數(shù)為

(4)

式中：Rr(t,xt,S)表示應(yīng)力S下突發(fā)失效的可靠度函數(shù)，xt表示突發(fā)失效發(fā)生時刻t對應(yīng)的退化量值，即X(t)=xt，τ為時間t的積分變量.

因此，競爭失效模式下t時刻產(chǎn)品的可靠度為

(5)

如果突發(fā)失效與退化量相互獨立，此時突發(fā)失效率可以表示為

λr(t,xt,S)=λr0(t,S)

(6)

式中：λr0(t)是不考慮性能退化影響情況下突發(fā)失效的失效率.

此時產(chǎn)品的可靠度為

(7)

以下研究僅考慮突發(fā)失效與退化量相關(guān)的情況.當(dāng)突發(fā)與退化不相關(guān)時，這時可采用傳統(tǒng)的構(gòu)造分布函數(shù)來刻畫突發(fā)失效的失效模型，再利用極大似然即可得到突發(fā)失效模型中的分布參數(shù).

根據(jù)競爭失效模型，在獲得應(yīng)力水平Si下競爭失效的可靠度函數(shù)后，可令應(yīng)力水平為S0，即可得到在常應(yīng)力S0下產(chǎn)品的各項可靠度估計.

1.3 加速退化方程加速退化方程[13]定義如下：記ζ(t)為產(chǎn)品退化量x(t)的某特征值或退化過程中的某參數(shù)，方程ζ(t)=A(t)eb(t)φ(S)或lnζ(t)=a(t)+b(t)φ(S)被稱為加速退化方程.其中S為應(yīng)力水平，φ(S)是應(yīng)力的函數(shù)，a(t)=ln(A(t))和b(t)是與退化量和應(yīng)力無關(guān)的參數(shù).

在實際工程應(yīng)用中，當(dāng)溫度作為加速應(yīng)力時，φ(S)=1/S.當(dāng)電應(yīng)力作為加速應(yīng)力時，φ(S)=lnS.

1.4 基于比例危險的突發(fā)失效模型為了對恒定應(yīng)力加速退化當(dāng)中突發(fā)失效進行更好的描述，本研究給出基于比例危險的突發(fā)失效模型.在恒定應(yīng)力加速退化試驗當(dāng)中，比例危險(proportional hazards，PH)族[11]是具有如下性質(zhì)的模型：同一時刻不同退化程度的兩個個體的失效率函數(shù)(失效率函數(shù)也被稱為危險函數(shù))與時間無關(guān)，僅與退化程度和當(dāng)時應(yīng)力有關(guān).設(shè)應(yīng)力Si下時刻t兩個產(chǎn)品的退化量分別為x1和x2，則這兩個產(chǎn)品在時刻t的失效率之比λ(t,x1,S)/λ(t,x2,S)與時間t無關(guān).對于突發(fā)失效時間T來說，如果其失效率函數(shù)具有如上性質(zhì)，則其條件失效率函數(shù)可寫成

λ(t,xt,S)=λ0(t)h(xt,S)

(8)

式中：λ0(t)是不考慮性能退化和應(yīng)力水平影響下突發(fā)失效的失效率，也叫做基準失效率函數(shù)；h(xt,S)是性能退化水平和應(yīng)力水平的函數(shù).

在該試驗下，可令h(xt,S)=exp(α1+α2xt+α3S)，式中：α1、α2、α3為未知回歸系數(shù).

對于上述基于比例危險的突發(fā)失效模型可以看到，獲得基準失效率函數(shù)即λ0(t)是該模型的關(guān)鍵所在，估計λ0(t)之后可以進一步得到突發(fā)失效的失效率函數(shù)，從而對競爭失效問題得到更清晰的認識.

2 競爭失效場合恒定應(yīng)力加速退化統(tǒng)計分析方法

2.2 突發(fā)與退化模型的參數(shù)估計

2.2.1 退化失效模型中的參數(shù)估計 在只考慮退化失效模式下，根據(jù)性能退化數(shù)據(jù)的散點圖特征，利用擬合優(yōu)度檢驗選取最合適的性能退化分布函數(shù).常用的退化量分布函數(shù)有對數(shù)正態(tài)分布、威布爾分布、正態(tài)分布等分布類型.

當(dāng)退化量分布為位置-刻度族模型時，一般地，位置參數(shù)可表示為退化時間(或其變換形式)、加速應(yīng)力(或其變換形式)的線性函數(shù)：

μ(t,S;β)=β1+β2t+β3φ(S)

(9)

當(dāng)退化量分布不為位置-刻度族模型時，則要通過加速退化方程確定退化過程某參數(shù)的函數(shù)形式.假設(shè)ω(t)為該分布的時變參數(shù)，由加速退化方程可得lnω(t)=a(t)+b(t)φ(S)，這里令a(t)=β1+β2t，b(t)=β3，則可將加速退化方程改寫為：

lnω(t)=β1+β2t+β3φ(S)

(10)

該加速退化方程同樣適用于威布爾分布，對數(shù)正態(tài)分布等位置-刻度族分布類型.

則樣本的聯(lián)合對數(shù)似然函數(shù)為：

(11)

對上述聯(lián)合對數(shù)似然函數(shù)采用極大似然估計可得到各個參數(shù)的估計值.

在獲得每個參數(shù)的估計后，則可得到產(chǎn)品退化失效模式下的可靠度Rd(t,xt,S).

2.2.2 比例危險模型λ0(t)的參數(shù)估計 當(dāng)突發(fā)失效與退化量相關(guān)時，引入基于比例危險的突發(fā)失效模型評估突發(fā)失效的失效率函數(shù).也即，在不同應(yīng)力下，對突發(fā)失效時間進行分布擬合，選擇最優(yōu)的分布類型作為基準分布類型，進一步確定基準失效率，從而表示出突發(fā)失效的失效率函數(shù).

首先，記錄產(chǎn)品在應(yīng)力Si下發(fā)生突發(fā)失效時對應(yīng)的退化量，記第l(l=1,2,…,Ni)個突發(fā)失效樣品失效時對應(yīng)的性能退化量為xl，該性能退化量對應(yīng)的突發(fā)失效時間為tl.

選擇合適的分布函數(shù)來擬合產(chǎn)品的突發(fā)失效時間，對于各個應(yīng)力上突發(fā)失效時間分布類型的確定，首先確定其備選分布類型.然后采用分布假設(shè)檢驗，從備選分布中確定各應(yīng)力下突發(fā)失效時間的最優(yōu)分布類型.

(12)

式中：Fn(x)為經(jīng)驗分布函數(shù)；F(x)為假設(shè)的累積分布函數(shù).

(13)

在確定各個應(yīng)力下的最優(yōu)分布類型后，選取各個應(yīng)力下的最優(yōu)分布類型為基準分布，對應(yīng)的失效率即為基準失效率.

2.2.3 比例危險模型h(xt,S)的參數(shù)估計 在選取基準失效率函數(shù)后，突發(fā)失效的失效率函數(shù)可表示為：

λ(t,xt,S)=λ0(t){exp(α1+α2xt+α3φ(S))}

(14)

式中：λ0(t)為基準失效率；α1、α2、α3為待估參數(shù)；xt為突發(fā)失效發(fā)生時刻t對應(yīng)的退化量.

根據(jù)失效率與可靠度的關(guān)系，可得到下列表達式：

(15)

當(dāng)突發(fā)失效時刻t已知時，其對應(yīng)的退化量xt隨即已知.觀察上式可知exp{α1+α2xt+α3φ(S)}與積分變量t無關(guān)，因此可將上式改寫為

(16)

將(16)式兩邊進行變換，可得到如下線性關(guān)系：

(17)

yt=α1+α2xt+α3φ(S)

(18)

對(18)式進行線性回歸即可得到α1、α2、α3的估計值.

而對(18)式進行線性回歸的前提是yt為已知量，但式中R(t,xt,S)卻是未知的，因此需要建立突發(fā)失效時間與可靠度的關(guān)系，下面將結(jié)合概率論中概率的基本概念以及古典概型等理論來建立產(chǎn)品突發(fā)失效時間與可靠度之間的聯(lián)系，進而得到y(tǒng)t與t的關(guān)系.

2.2.4 突發(fā)失效可靠度函數(shù)預(yù)計 為了簡化公式，以下均在應(yīng)力水平Si下計算突發(fā)失效故障率并進行參數(shù)估計.

首先，從退化數(shù)據(jù)中判斷發(fā)生突發(fā)失效的產(chǎn)品并記錄其發(fā)生突發(fā)失效時對應(yīng)的突發(fā)失效時間.記第l(l=1,2,…,Ni)個突發(fā)失效樣品失效時對應(yīng)的突發(fā)失效時間為tl，并將其進行升序排列，得到{tn}，其中n=1,2,…,Ni.

建立突發(fā)失效模式下突發(fā)失效故障率的函數(shù)關(guān)系，即當(dāng)突發(fā)失效時間為tn時，突發(fā)失效故障率為：

(19)

由突發(fā)失效故障率與可靠度的關(guān)系，便可得到突發(fā)失效的可靠度.

綜上，可得到基于比例危險的突發(fā)失效模型的失效率函數(shù)，進而結(jié)合可得到競爭失效下恒定應(yīng)力加速退化產(chǎn)品的可靠度函數(shù).

3 算例

現(xiàn)有一批激光器進行退化試驗，該試驗利用了GaAs激光器的工作電流，在80 ℃溫度下隨時間變化的百分比數(shù)據(jù)記錄[14].共有22個樣本參加試驗，每隔500 h測試一次樣本，至4 000 h為止.當(dāng)激光器的工作電流增加至其額定電流的110%時，則認為該激光器發(fā)生退化失效，每個激光器監(jiān)測8次.其中15個樣本發(fā)生性能退化，另外7個樣本發(fā)生突發(fā)失效.現(xiàn)利用本文中提出的方法對該型激光器進行可靠性評估.并與文獻[8]中所給的可靠度曲線進行對比，驗證本研究所建模型的正確性.

1)退化失效模型中的參數(shù)估計.

由文獻[14]可知，1、2、4、6、7、8、10、11、13、14、16、18、19、20、22號樣本為性能退化數(shù)據(jù).圖1顯示的是該型激光器性能退化數(shù)據(jù)的變化趨勢.

圖1 激光器性能退化數(shù)據(jù)

通過分析可知各測量時刻的退化量服從威布爾分布故采用威布爾分布進行退化失效建模.由加速退化方程(10)式可得到恒定應(yīng)力下威布爾分布函數(shù)為

(20)

式中：形狀參數(shù)m為常數(shù)，尺度參數(shù)ln[η(t,S)]=β1+β2t+β3φ(S)，β1、β2、β3為未知參數(shù).

鑒于該激光器試驗僅在額定應(yīng)力80 ℃下進行，所以可將上述尺度參數(shù)進行變換以滿足試驗，但并不影響該方法在恒定應(yīng)力加速退化試驗當(dāng)中應(yīng)用.由于溫度應(yīng)力S為常數(shù)，可將β1與β3φ(S)合并，令β0=β1+(1/80)β3，尺度參數(shù)則變形為ln[η(t)]=β0+β2t.因此，激光器退化量的分布函數(shù)為

(21)

則所有樣本的對數(shù)似然函數(shù)可表示為

(22)

(23)

2)比例危險模型λ0(t)的參數(shù)估計.

由文獻[14]可知，3、5、9、12、15、17、21號樣本為突發(fā)失效數(shù)據(jù)，突發(fā)失效時間及其對應(yīng)的退化量見表1.

表1 突發(fā)失效時間及其退化量

確定突發(fā)失效時間分布模型，備選分布類型有：正態(tài)分布、Weibull分布、對數(shù)正態(tài)分布、Gamma分布等.鑒于鑒于該激光器試驗僅在額定應(yīng)力80 ℃下進行，所以在應(yīng)力80 ℃下采用AD檢驗法確定的最優(yōu)分布即為基準分布.計算80 ℃的檢驗統(tǒng)計量AD值(見表2).

表2 80 ℃下檢驗統(tǒng)計量AD值

(24)

式中：φ(·)與Φ(·)分別指標(biāo)準正態(tài)分布的密度函數(shù)和分布函數(shù).

3)突發(fā)失效可靠度函數(shù)預(yù)計與比例危險模型h(xt,S)的參數(shù)估計.

將表1中的突發(fā)失效數(shù)據(jù)按突發(fā)失效時間排序，并按照1.2.4求出突發(fā)失效時間對應(yīng)的突發(fā)失效故障率、可靠度，如表3.

表3 80 ℃下突發(fā)失效相關(guān)數(shù)據(jù)

鑒于該激光器試驗僅在額定應(yīng)力80 ℃下進行，所以可將(18)式進行變換以滿足試驗，但并不影響該方法在恒定應(yīng)力加速退化試驗中的應(yīng)用.由于溫度應(yīng)力S為常數(shù)，可將α1與α3φ(S)合并，令α0=α1+(1/80)α3，則(18)式變形為

yt=α0+α2xt

(25)

由此得到基于比例危險的突發(fā)失效模型的失效率函數(shù)為

(26)

依據(jù)提出的可靠性建模方法，得到競爭失效情況下的可靠度為

(27)

圖2 競爭失效可靠度曲線

由圖2可見：在監(jiān)測初期，激光器突發(fā)失數(shù)據(jù)較少，M1所得曲線略高于M2，但總的來說，兩條可靠度曲線相差不大，與實際退化情況基本相符.在監(jiān)測后期，突發(fā)失效數(shù)據(jù)增多，并且退化失效數(shù)據(jù)也開始增多.以4 000 h為例，文獻[8]中給出的可靠度為0.65，本研究給出的可靠度為0.46，而通過頻率去估計該時間點的可靠度為0.54，相較于文獻[8]，本研究的可靠度值較為接近，更能客觀地反映這一現(xiàn)象.相比M2，M1得到的可靠性評估精度略優(yōu)且略保守，符合工程上的做法.除此之外，采用基于比例危險的突發(fā)失效模型，可以更方便地推廣到競爭失效下恒定應(yīng)力場合，適用性較文獻[8]中的結(jié)果更好.

4 結(jié)論

為了解決恒定應(yīng)力加速退化試驗中競爭失效可靠性評估問題，采用加速退化方程和比例危險模型對該問題進行刻畫，豐富了該方面的研究.

本研究基于突發(fā)與退化競爭失效模型，對突發(fā)失效和退化失效分別展開研究.

1)針對偽失效壽命方法的不足，采用退化量分布函數(shù)描述產(chǎn)品性能的退化過程.

2)對產(chǎn)品在性能退化中伴有的突發(fā)失效，采用基于比例危險的突發(fā)失效模型，通過AD檢驗法獲得基準失效率函數(shù)λ0(t)，而對比例危險模型中的h(xt,S)，則通過非參數(shù)方法和線性回歸得到h(xt,S)中的參數(shù)估計.

本文中提出的競爭失效場合下恒定應(yīng)力加速退化可靠性模型，更加客觀地反映了退化與突發(fā)的競爭失效過程，具有一定的工程應(yīng)用價值，但對于含有多種退化失效和突發(fā)失效的競爭失效問題還需要繼續(xù)研究.