魯 相，汪亞順，陳 循，譚源源

(1. 國防科技大學(xué) 智能科學(xué)學(xué)院， 湖南 長沙 410073; 2. 中國空氣動力研究與發(fā)展中心 設(shè)備設(shè)計(jì)及測試技術(shù)研究所， 四川 綿陽 621000)

隨著產(chǎn)品質(zhì)量水平的不斷提高，近代研制出的產(chǎn)品普遍具有高可靠、長壽命的特點(diǎn)，例如發(fā)光二極管及橡膠密封件等元器件中位壽命可以超過10年[1]。對于這類高可靠性產(chǎn)品，現(xiàn)代工業(yè)廣泛利用加速退化試驗(yàn)預(yù)測其使用壽命和貯存壽命。在加速退化試驗(yàn)中，將產(chǎn)品置于加速應(yīng)力水平下并記錄產(chǎn)品樣本的性能退化數(shù)據(jù)，利用加速模型建立應(yīng)力與退化速率之間的關(guān)系。當(dāng)各試驗(yàn)應(yīng)力水平下的失效機(jī)理相同時，就可以將加速應(yīng)力下的試驗(yàn)結(jié)果外推，進(jìn)而估計(jì)在正常應(yīng)力狀況下的產(chǎn)品壽命。各試驗(yàn)應(yīng)力水平下失效機(jī)理具有一致性是加速退化試驗(yàn)結(jié)論具備有效性的首要前提，因此需要研究加速退化試驗(yàn)機(jī)理一致性判別方法。

利用試驗(yàn)數(shù)據(jù)統(tǒng)計(jì)分析方法或物理分析測試手段，判斷失效機(jī)理在試驗(yàn)各加速應(yīng)力水平下是否發(fā)生變化，這個過程被定義為失效機(jī)理的一致性判別。目前失效機(jī)理一致性判別方法可分為3類：概率圖檢驗(yàn)法[2-4]、基于模型參數(shù)一致性檢驗(yàn)法[5-9]、基于灰色預(yù)測模型檢驗(yàn)法[10-11]。這些一致性判別方法各有優(yōu)劣，但存在共同的缺點(diǎn)：通常只對壽命服從某一特定分布的產(chǎn)品有效，應(yīng)用場合受限，因而難以推廣到其他統(tǒng)計(jì)分布，并且這些方法對判別結(jié)果的誤判風(fēng)險沒有開展定量分析。因此，提出基于似然比檢驗(yàn)的失效機(jī)理一致性判別方法，該方法適用于多個統(tǒng)計(jì)分布，可以定量分析誤判風(fēng)險，并且模型參數(shù)估計(jì)精度較高。

1 機(jī)理一致性判別的數(shù)學(xué)模型及參數(shù)估計(jì)

1.1 兩類加速模型

通常利用混合效應(yīng)模型[12-13]對高可靠產(chǎn)品退化建模，并認(rèn)為性能退化量初始值B及退化軌跡形狀參數(shù)α是固定不變的，不同個體間的退化速率存在差異。產(chǎn)品的退化速率Kij在各應(yīng)力水平Si下一般服從對數(shù)正態(tài)分布LN(μi,σ2),其中σ為對數(shù)標(biāo)準(zhǔn)差，μi為對數(shù)均值(與應(yīng)力水平有關(guān))。選擇如式(1)所示的混合效應(yīng)模型。

(1)

在混合效應(yīng)模型下，當(dāng)失效機(jī)理不變時，提高載荷或溫度的試驗(yàn)產(chǎn)品統(tǒng)一加速模型為：

μi=a+b·φ(Si)

(2)

式中，φ(Si)為轉(zhuǎn)換應(yīng)力水平，溫度為加速應(yīng)力時φ(Si)=1/Si，載荷為加速應(yīng)力時φ(Si)=lnSi。

當(dāng)加速應(yīng)力范圍較大導(dǎo)致失效機(jī)理變化時，式(2)中的參數(shù)a和b不再為常數(shù)，式(2)在較高應(yīng)力水平下無效。比式(2)更適用的模型可以表述為：

μi=ai+bi·φ(Si)

(3)

可得到如下原假設(shè)H0:μi=a+b·φ(Si)和備擇假設(shè)H1:μi=ai+bi·φ(Si)。

從上述定義可以看出，在原假設(shè)中假定退化速率-應(yīng)力關(guān)系為對數(shù)線性，將原假設(shè)對應(yīng)模型稱為對數(shù)線性加速模型。相比之下，在備擇假設(shè)中模型參數(shù)與應(yīng)力水平有關(guān)，將備擇假設(shè)對應(yīng)模型稱為非對數(shù)線性加速模型。式(2)與式(3)的表達(dá)式表明原假設(shè)對應(yīng)模型是備擇假設(shè)對應(yīng)模型的一種特例。在對數(shù)線性加速模型中參數(shù)a和b被限定為常數(shù)，導(dǎo)致其獨(dú)立參數(shù)少于非對數(shù)線性加速模型。如果兩種假設(shè)模型中一種模型是另一模型的特例，則可用似然比檢驗(yàn)比較兩種模型的擬合優(yōu)度[14-15]，因此似然比檢驗(yàn)可用于檢驗(yàn)上述原假設(shè)。

1.2 似然比檢驗(yàn)的基本準(zhǔn)則

假定有來自于分布T的隨機(jī)樣本t1,t2,…,tn，令θ表示該分布未知參數(shù)向量，則該樣本的似然函數(shù)為：

(4)

令L0(θ)為對數(shù)線性加速模型參數(shù)的似然函數(shù)，L1(θ)為非對數(shù)線性加速模型參數(shù)的似然函數(shù)。對數(shù)線性加速模型參數(shù)的極大似然估計(jì)(Maximum Likelihood Estimation, MLE)為：

(5)

非對數(shù)線性加速模型參數(shù)的極大似然估計(jì)為：

(6)

似然比統(tǒng)計(jì)量可表示為：

(7)

由于包含更多獨(dú)立參數(shù)，通常非對數(shù)線性模型對樣本觀測值的擬合效果比對數(shù)線性模型好，但兩類模型也可能都有較好的擬合效果。非對數(shù)線性模型參數(shù)的似然函數(shù)值不小于對數(shù)線性模型的似然函數(shù)值，似然比統(tǒng)計(jì)量滿足：

(8)

如果對數(shù)線性模型與非對數(shù)線性模型一樣適用，則似然比較大，反之較小。當(dāng)似然比統(tǒng)計(jì)量數(shù)值過小時似然比檢驗(yàn)法會拒絕原假設(shè)(對數(shù)線性加速模型)，似然比檢驗(yàn)法的臨界域或拒絕域?yàn)椋?/p>

W={λ≤c} 0

(9)

式中，臨界值c與參數(shù)λ、β有關(guān)。β是似然比檢驗(yàn)中規(guī)定的顯著性水平,β表示將失效機(jī)理相同誤判為失效機(jī)理變化的風(fēng)險大小。

1.3 數(shù)學(xué)模型及參數(shù)估計(jì)

根據(jù)上述分析，首次應(yīng)用似然比檢驗(yàn)法辨識加速退化試驗(yàn)中加速模型的變化，并計(jì)算兩種假設(shè)模型下對數(shù)似然函數(shù)的極大值。設(shè)恒定應(yīng)力加速退化試驗(yàn)(Constant Stress Accelerated Degradation Testing, CSADT)采用q個應(yīng)力水平,在每個加速應(yīng)力水平Si(i=1,2,…,q)下測試ni個試驗(yàn)樣本的性能退化量，測試mi次，測試的時間節(jié)點(diǎn)為tik(k=1,2,…,mi)，在每個應(yīng)力水平下的試驗(yàn)截止時間為ti，測得退化量數(shù)據(jù)為yijk(j=1,2,…,ni)。建模時取樣本的性能退化量對數(shù)值lnyijk統(tǒng)計(jì)分析，混合效應(yīng)模型下其條件累積分布及條件概率密度函數(shù)可以表示為：

(10)

(11)

其中，Φ(·)和φ(·)表示標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)及概率密度函數(shù)。由式(10)及式(11)可得恒定應(yīng)力加速退化試驗(yàn)數(shù)據(jù)的似然函數(shù)為：

(12)

式中，f(Kij)代表退化速率Kij的概率密度函數(shù)。

(13)

假定采用溫度作為加速應(yīng)力，在對數(shù)線性加速模型中ai=a,bi=b，取向量δ=(lnB,α,a,b)，Kij=exp(a+b/Si+cij)，yijk=(yij1,…,yijmi)T，tik=(ti1,…,timi)T，協(xié)方差矩陣D=σ2，cij服從正態(tài)分布N(0,σ2)。采用近似極大似然估計(jì)法求解模型參數(shù)。

步驟1：在固定協(xié)方差矩陣的情況下，用非線性最小二乘法獲得加速模型參數(shù)a、b、cij和形狀參數(shù)α及退化量初始值B的估計(jì)值，估計(jì)值是使式(14)最小化的解。

(14)

步驟2：首先將模型函數(shù)在第一步獲得的各參數(shù)估計(jì)值處作Taylor展開，根據(jù)模型函數(shù)表達(dá)式，有

其中， *代表哈達(dá)馬積。

然后利用第一步獲得的估計(jì)值求取協(xié)方差矩陣D及其他參數(shù)的極大似然函數(shù)表達(dá)式，取

可得對數(shù)線性加速模型的對數(shù)似然函數(shù)為：

(15)

設(shè)定好數(shù)值迭代算法的初始解，利用Newton-Raphson算法[16-19]求解未知參數(shù)估計(jì)值，按式(15)計(jì)算對數(shù)似然函數(shù)的極大值。

在非對數(shù)線性加速模型中，取向量κ=(lnB,α,μ1,…,μq)，Kij=exp(μi+cij)，yijk=(yij1,…,yijmi)T，tik=(ti1,…,timi)T，協(xié)方差矩陣D=[σ2]。

按以下步驟進(jìn)行參數(shù)估計(jì)。

步驟1：在固定協(xié)方差矩陣的情況下，用非線性最小二乘法獲得加速模型參數(shù)μi、cij和形狀參數(shù)α及退化量初始值B的估計(jì)值，估計(jì)值是使式(16)最小化的解。

(16)

步驟2：首先將模型函數(shù)在第一步獲得的各參數(shù)估計(jì)值處作Taylor展開，根據(jù)模型函數(shù)表達(dá)式，有

同樣求取協(xié)方差矩陣D及其他參數(shù)的極大似然函數(shù)表達(dá)式，取

可得非對數(shù)線性加速模型的對數(shù)似然函數(shù)為：

(17)

利用Newton-Raphson算法類似得到對數(shù)似然函數(shù)的極大值。

2 失效機(jī)理一致性判別準(zhǔn)則及方法流程

算出兩種假設(shè)下似然函數(shù)極大值后，需要確定似然比統(tǒng)計(jì)量的分布，以構(gòu)建判別準(zhǔn)則。多數(shù)情況下難以確定似然比統(tǒng)計(jì)量λ的精確抽樣分布，但可以確定大樣本場合下對數(shù)似然比統(tǒng)計(jì)量的漸進(jìn)分布，令

(18)

根據(jù)Wilks提出的廣義似然比統(tǒng)計(jì)量的極限分布定理[20]，在原假設(shè)成立的情況下，統(tǒng)計(jì)量Λ漸進(jìn)服從自由度為v的卡方分布χ2(v)。自由度v等于在假設(shè)H0和H1下的獨(dú)立參數(shù)數(shù)目之差，有v=q-2。由式(9)可知，當(dāng)λ≤c時原假設(shè)H0被拒絕，因此當(dāng)Λ≥-2lnc時對數(shù)線性加速模型無效。利用顯著性水平β及統(tǒng)計(jì)量Λ的分布，可以得出

(19)

式中，等式右端表示卡方分布χ2(q-2)的1-β分位數(shù)。

利用式(15)和式(17)計(jì)算恒定應(yīng)力加速退化試驗(yàn)中兩類模型的極大對數(shù)似然函數(shù)的差值，得到：

ΛCSM=-2[lnL0(δ,D,σε)-lnL1(κ,D,σε)]

(20)

失效機(jī)理一致性的判別準(zhǔn)則是：

1)如果0<ΛCSM<-2lnc，失效機(jī)理在所有應(yīng)力水平下相同；

2)如果ΛCSM≥-2lnc，失效機(jī)理在某一高應(yīng)力水平下變化。

根據(jù)上述分析，基于似然比檢驗(yàn)的失效機(jī)理一致性判別方法流程如圖1所示。

3 仿真算例

在發(fā)光二極管(Light Emitting Diode, LED)加速試驗(yàn)中，最高溫度下激活能可能發(fā)生變化。本例中設(shè)定退化失效機(jī)理發(fā)生變化，在378 K激活能從16.1 kJ/mol轉(zhuǎn)變?yōu)?2.2 kJ/mol。

圖1 失效機(jī)理一致性判別方法流程Fig.1 Flow chart of identification method of failure mechanism consistency

仿真試驗(yàn)溫度應(yīng)力水平為298 K、338 K、378 K，試驗(yàn)樣本量為6，測量間隔為336 h，每個水平下測量次數(shù)為11。設(shè)定退化模型參數(shù)分別為a1=-0.546 6，b1=-1933.25(298～338 K)，a2=5.173 1，b2=-3866.49(338～378 K)，B=1.051 6，α=0.65，σ=0.185，σε=0.014。利用式(1)生成仿真退化軌跡。

利用近似極大似然估計(jì)法分析仿真數(shù)據(jù)，得出對數(shù)線性模型中未知參數(shù)的極大似然估計(jì)，將極大似然估計(jì)值與未知參數(shù)真值進(jìn)行比較，結(jié)果如表1所示。

表1結(jié)果表明，采用指定的加速模型進(jìn)行建模及統(tǒng)計(jì)分析時，加速模型參數(shù)估計(jì)值a、b與真值的偏差較大，估計(jì)值近似于兩真值的均值，導(dǎo)致在298 K以下的預(yù)測壽命偏大。

表1 對數(shù)線性模型參數(shù)估計(jì)值

非對數(shù)線性模型中未知參數(shù)極大似然估計(jì)值與真值的比較如表2所示，可以看出，模型各參數(shù)的估計(jì)值與真值相差較小，非對數(shù)線性模型擬合效果明顯優(yōu)于對數(shù)線性模型。

進(jìn)一步可以算得溫度水平在298～378 K下檢驗(yàn)統(tǒng)計(jì)量為：

ΛCSM=-2lnλ=-2(lnL0-lnL1)

=-2×(510.254 7-518.963 2)=17.417

(21)

在298～378 K溫度范圍內(nèi)的試驗(yàn)應(yīng)力水平數(shù)為3，給定顯著性水平β=5%，臨界值滿足-2lnc=3.84<ΛCSM。因此在298～378 K范圍內(nèi)失效機(jī)理發(fā)生變化，這一判別結(jié)果與失效機(jī)理發(fā)生變化的原始設(shè)定相符，驗(yàn)證了方法的正確性。

表2 非對數(shù)線性模型參數(shù)估計(jì)值

4 應(yīng)用實(shí)例

丁腈橡膠密封件是典型的高可靠性產(chǎn)品，但易于發(fā)生熱氧老化。受熱氧老化的影響，丁腈橡膠密封件在貯存或使用期密封性能會不斷退化。隨著密封性能的退化，橡膠密封性能的高彈性復(fù)原力不斷減小，壓縮永久變形不斷增大。因此可以通過觀測壓縮永久變形cs的變化來監(jiān)測退化趨勢，當(dāng)壓縮永久變形上升到某一臨界值時，認(rèn)為橡膠密封件密封失效。試驗(yàn)采用安裝于系統(tǒng)中的橡膠密封件的壓縮率，樣本可以代表在空氣介質(zhì)存儲下系統(tǒng)內(nèi)受壓丁腈橡膠密封件的退化。試驗(yàn)總樣本量為20，選擇在5個溫度應(yīng)力水平(333 K、343 K、353 K、373 K、393 K)下進(jìn)行加速退化試驗(yàn)，每個應(yīng)力水平下樣本量為4。密封件在333 K及393 K應(yīng)力水平下測得的壓縮永久變形量如圖2所示，其他應(yīng)力水平下曲線形狀類似，不再給出。根據(jù)圖1所示流程對數(shù)據(jù)進(jìn)行失效機(jī)理一致性判別。建立退化模型，將測得的壓縮永久變形cs(無量綱)轉(zhuǎn)換為性能退化數(shù)據(jù)y，轉(zhuǎn)換關(guān)系式為y=1-cs，采用混合效應(yīng)模型(1)對試樣退化規(guī)律進(jìn)行描述。

(a) 333 K丁腈密封圈加速退化試驗(yàn)數(shù)據(jù) (a) CSADT data of nitrile seal rings at 333 K

(b) 393 K丁腈密封圈加速退化試驗(yàn)數(shù)據(jù) (b) CSADT data of nitrile seal rings at 393 K圖2 丁腈密封圈恒定應(yīng)力加速退化試驗(yàn)數(shù)據(jù)Fig.2 CSADT data of nitrile seal rings

利用近似極大似然估計(jì)方法，令ni=4(i=1,2,3)，可求解333～353 K、333～373 K及333～393 K下模型的未知參數(shù)，設(shè)定顯著性水平β=10%，進(jìn)一步算得密封圈試樣在三個不同溫度范圍的假設(shè)檢驗(yàn)結(jié)果，如表3所示。

表3 密封圈失效機(jī)理一致性判別

根據(jù)失效機(jī)理一致性判別方法，結(jié)合表3計(jì)算結(jié)果可知，試驗(yàn)所用丁腈橡膠密封圈在373 K以上時溫度失效機(jī)理發(fā)生變化。Wise等[21]、Le等[22]對類似產(chǎn)品開展失效物理分析，得到368 K溫度以上丁腈橡膠失效機(jī)理變化的結(jié)果，與本文所得結(jié)論基本一致，驗(yàn)證了本文方法的有效性。

5 結(jié)論

高可靠性產(chǎn)品在進(jìn)行加速退化試驗(yàn)時，失效機(jī)理在高應(yīng)力水平下可能發(fā)生變化，導(dǎo)致錯誤的壽命預(yù)測結(jié)果?；谒迫槐葯z驗(yàn)的加速退化試驗(yàn)機(jī)理一致性判別方法，適用范圍廣，能夠定量分析誤判風(fēng)險，并給出機(jī)理一致性判別結(jié)論，得到準(zhǔn)確的壽命預(yù)測結(jié)果。建立非線性加速模型與高可靠長壽命產(chǎn)品失效機(jī)理變化間的聯(lián)系，并將似然比檢驗(yàn)應(yīng)用于加速退化試驗(yàn)中對加速模型的辨識，建立了有效的失效機(jī)理一致性判別準(zhǔn)則及方法。

通過仿真案例及應(yīng)用實(shí)例，進(jìn)一步驗(yàn)證了所提方法的有效性。多應(yīng)力加速退化試驗(yàn)下產(chǎn)品失效機(jī)理的一致性判別較為復(fù)雜，相關(guān)方法還需要進(jìn)一步研究。