基于隨機(jī)相關(guān)的電子部件二元加速退化可靠性評(píng)估

蓋炳良， 滕克難， 王浩偉， 王文雙， 陳健， 宦婧

(1. 海軍航空大學(xué) 煙臺(tái) 264001； 2. 中國(guó)人民解放軍91115部隊(duì), 舟山 316000； 3. 江蘇科技大學(xué)計(jì)算機(jī)學(xué)院, 鎮(zhèn)江 212003)

隨著材料科學(xué)與生產(chǎn)制造工藝的發(fā)展進(jìn)步，軍工領(lǐng)域出現(xiàn)了越來(lái)越多的高可靠性長(zhǎng)壽命電子產(chǎn)品。此類(lèi)產(chǎn)品在額定工作應(yīng)力下性能退化緩慢，僅僅依靠性能檢測(cè)數(shù)據(jù)已經(jīng)難以快速有效地進(jìn)行可靠性評(píng)估。因而加速退化試驗(yàn)成為可靠性評(píng)估的有效途徑。通常認(rèn)為溫度是導(dǎo)致電子部件退化失效的最主要環(huán)境應(yīng)力類(lèi)型，溫度對(duì)多種失效機(jī)理具有較好的加速效果，如電子元器件的二次俘獲、表面電荷擴(kuò)散、電遷移等，因而常將溫度作為加速應(yīng)力。

通常認(rèn)為，由于生產(chǎn)工藝、貯存條件、人員操作等多方面主客觀(guān)因素，產(chǎn)品的個(gè)體差異性是客觀(guān)存在的。目前，考慮個(gè)體差異的隨機(jī)過(guò)程模型在可靠性領(lǐng)域得到了較多研究應(yīng)用，如最常用的3類(lèi)隨機(jī)過(guò)程模型：Wiener過(guò)程[1-2]、Gamma過(guò)程[3]、Inverse Gaussian過(guò)程模型[4-5]等，但以上研究都是在額定應(yīng)力下的產(chǎn)品可靠性評(píng)估。如若考慮溫度加速應(yīng)力作用下的電子部件二元退化建模，考慮個(gè)體差異的隨機(jī)過(guò)程可靠性評(píng)估則需要解決以下2個(gè)重要問(wèn)題：一是如何確定此類(lèi)模型參數(shù)與加速應(yīng)力的約束關(guān)系；二是如何進(jìn)行相關(guān)性建模。

確定模型參數(shù)與加速應(yīng)力的約束關(guān)系是基于加速退化數(shù)據(jù)可靠性評(píng)估的前提，也是一個(gè)研究難點(diǎn)。目前，采用的模型通常都是隨機(jī)過(guò)程基本模型，即針對(duì)的都是不考慮個(gè)體差異的隨機(jī)過(guò)程模型參數(shù)與加速應(yīng)力的約束關(guān)系。此種情況下，解決方法大多是根據(jù)經(jīng)驗(yàn)假定模型參數(shù)與加速應(yīng)力的關(guān)系，比如Wiener過(guò)程模型，就存在著2個(gè)不同假定：假定漂移參數(shù)與應(yīng)力有關(guān)而擴(kuò)散參數(shù)與加速應(yīng)力無(wú)關(guān)[6-9]；假定漂移參數(shù)和擴(kuò)散參數(shù)都與加速應(yīng)力相關(guān)[10-11]等。顯然，這易導(dǎo)致評(píng)估結(jié)果可信度不高。另一種解決方法就是通過(guò)加速因子不變?cè)瓌t推導(dǎo)模型參數(shù)與加速應(yīng)力的約束關(guān)系，加速因子不變?cè)瓌t是指為保證加速因子具有工程應(yīng)用性，則加速因子應(yīng)當(dāng)是一個(gè)不隨時(shí)間變化的常數(shù)，比如文獻(xiàn)[12-13]推導(dǎo)得到了Wiener過(guò)程模型漂移參數(shù)和擴(kuò)散參數(shù)與加速應(yīng)力的約束關(guān)系。以上，都是加速應(yīng)力作用下，針對(duì)不考慮個(gè)體差異的隨機(jī)過(guò)程模型的，而文獻(xiàn)[14]采用加速因子不變?cè)瓌t推導(dǎo)得到了加速應(yīng)力下，考慮個(gè)體差異的Wiener過(guò)程模型參數(shù)與應(yīng)力的約束關(guān)系。

二元加速退化建模中，相關(guān)性建模通常采用基于Copula函數(shù)的建模方法[15]，也有采用基于馬氏距離的二元降維建模方法[16]的。在基于Copula函數(shù)的建模方法中，基本都將Copula函數(shù)參數(shù)視為固定值，即忽略了產(chǎn)品個(gè)體之間存在的隨機(jī)相關(guān)性。文獻(xiàn)[17]提出了隨機(jī)相關(guān)的建模思想，但只應(yīng)用在額定應(yīng)力下，沒(méi)有進(jìn)一步考慮加速應(yīng)力作用下Copula函數(shù)參數(shù)與加速應(yīng)力的關(guān)系。

綜上，為解決電子部件二元加速退化可靠性評(píng)估難題，建立基于貝葉斯方法的二元隨機(jī)相關(guān)性能退化可靠性評(píng)估框架。采用Wiener RDV(Random Drift-Volatility)模型進(jìn)行邊緣分布建模，在此基礎(chǔ)上，構(gòu)建基于Copula函數(shù)的多個(gè)隨機(jī)相關(guān)模型，采用兩階段貝葉斯參數(shù)估計(jì)方法進(jìn)行參數(shù)估計(jì)(通過(guò)OpenBUGS軟件實(shí)現(xiàn)基于MCMC(Markov Chain Monte Carlo)的參數(shù)估計(jì))，進(jìn)行相關(guān)性模型選擇，最后基于蒙特卡羅方法進(jìn)行可靠度計(jì)算，并應(yīng)用實(shí)例驗(yàn)證了所提方法有效性。

1 二元隨機(jī)相關(guān)可靠性分析框架

基于Wiener過(guò)程的二元加速性能退化可靠性分析，以邊緣退化建模為基礎(chǔ)，利用Copula函數(shù)實(shí)現(xiàn)二元相關(guān)關(guān)系建模。主要包括3部分內(nèi)容：①邊緣分布建模和相關(guān)性建模；②兩階段貝葉斯參數(shù)估計(jì)；③基于蒙特卡羅方法的可靠性評(píng)估。首先，利用樣本性能退化數(shù)據(jù)，采用RDV模型建立邊緣分布模型，并考慮產(chǎn)品間具有隨機(jī)相關(guān)性，建立基于Copula函數(shù)的相關(guān)性模型；然后，采用兩階段貝葉斯參數(shù)估計(jì)方法對(duì)邊緣分布模型進(jìn)行參數(shù)估計(jì)，在此基礎(chǔ)上得到邊緣分布退化增量的累積概率分布值，將其作為Copula函數(shù)的輸入值，并通過(guò)Copula函數(shù)模型選擇，確定相對(duì)最優(yōu)的相關(guān)模型；最后，在參數(shù)估計(jì)和模型選擇基礎(chǔ)上，采用蒙特卡羅方法建立可靠度函數(shù)曲線(xiàn)，完成產(chǎn)品的可靠性評(píng)估?；究蚣苋鐖D1所示。

圖1 二元相關(guān)退化建?？蚣蹻ig.1 Framework of bivariate correlation degradation modeling

本文采用如下3個(gè)假定。

假定1同型產(chǎn)品，不僅在性能退化階段具有個(gè)體差異，而且在產(chǎn)品之間的相關(guān)關(guān)系上也可能具有個(gè)體差異性。

假定2同型產(chǎn)品相關(guān)關(guān)系的差異性，是指Copula函數(shù)參數(shù)的差異性，即同型產(chǎn)品具有相同類(lèi)型的Copula函數(shù)，但Copula函數(shù)的參數(shù)值不同。

假定3加速應(yīng)力可能會(huì)影響相關(guān)關(guān)系，但是應(yīng)力不改變Copula函數(shù)類(lèi)型，只改變Copula函數(shù)參數(shù)值。

2 Wiener RDV模型

文獻(xiàn)[14]給出了一種考慮個(gè)體差異的Wiener RDV模型，該模型將性能退化過(guò)程的均值和方差都確定為隨機(jī)變量，該模型描述為

(1)

式中：X(t)為退化量;B(·)為標(biāo)準(zhǔn)Brown運(yùn)動(dòng)函數(shù);μ為漂移參數(shù);σ為擴(kuò)散參數(shù);Λ(t)=tc為時(shí)間函數(shù)；a、b、e、h為超參數(shù)；N(h,e/?)表示均值為h，方差為e/?的正態(tài)分布。模型的累積分布函數(shù)為

(2)

其中：D為閾值；F2a(·)為自由度為2a的t分布函數(shù)。文獻(xiàn)[14]采用加速因子不變?cè)瓌t推導(dǎo)得到了該模型參數(shù)與加速應(yīng)力的約束關(guān)系：

(3)

式中：ak、bk、ck、ek和hk為應(yīng)力sk的模型參數(shù)；al、bl、cl、el和hl為應(yīng)力sl的模型參數(shù)；AF為應(yīng)力sk相對(duì)于應(yīng)力sl的加速因子。由式(3)可采用Arrhenius模型將溫度應(yīng)力s下參數(shù)表示為

(4)

式中：γRDV(1)、γRDV(2)、γRDV(3)和γRDV(4)為待定參數(shù)。

3 二元隨機(jī)相關(guān)加速退化數(shù)據(jù)建模

假定產(chǎn)品具有2個(gè)壽命表征參數(shù)X1和X2，邊緣性能退化過(guò)程都服從Wiener RDV模型。tijk為第i個(gè)樣本第k個(gè)參數(shù)第j次測(cè)量時(shí)刻，Xijk為tijk時(shí)刻測(cè)量值，k=1,2；i=1,2,…,N；j=1,2,…,M；N為樣本數(shù)量，M為測(cè)量次數(shù)。對(duì)所有樣本進(jìn)行了相同次數(shù)的測(cè)量，不同性能參數(shù)在相同時(shí)刻進(jìn)行測(cè)量。則易得性能退化數(shù)據(jù)結(jié)構(gòu)為

(5)

1) 模型A

不考慮相關(guān)關(guān)系個(gè)體差異性，Copula函數(shù)參數(shù)θ是一個(gè)固定參數(shù)，這是目前通常采用的模型，即

(6)

模型B～模型F的性能退化模型部分與模型A相同，因而以下內(nèi)容只描述相關(guān)性模型部分。

2) 模型B

(7)

3) 模型C

考慮相關(guān)關(guān)系存在個(gè)體差異性，Copula函數(shù)參數(shù)是服從正態(tài)分布的隨機(jī)變量，即

(8)

4) 模型D

(9)

5) 模型E

(10)

6) 模型F

(11)

4 參數(shù)估計(jì)及模型選擇

本文采用常用的Gaussian、Frank、Gumbel 和Clayton Copula函數(shù)進(jìn)行模型選擇，其分布函數(shù)及參數(shù)如表1所示。這4類(lèi)Copula函數(shù)描述了不同的相關(guān)關(guān)系結(jié)構(gòu)，Gaussian Copula和Frank Copula描述的是對(duì)稱(chēng)相關(guān)結(jié)構(gòu)，Gumbel Copula函數(shù)描述較強(qiáng)上尾相關(guān)特征的變量間的相關(guān)關(guān)系，Clayton Copula函數(shù)則描述具有較強(qiáng)下尾相關(guān)特征的變量相關(guān)關(guān)系。

Copula函數(shù)模型選擇方法，可以采用散點(diǎn)圖，偏差信息準(zhǔn)則(DIC)值以及Kendallτ的非參數(shù)估計(jì)值等方法。

表1 Copula 函數(shù)Table 1 Copula function

(12)

DIC值是常用的用于模型選擇的一種定量方法，定義為

(13)

5 基于蒙特卡羅的可靠度計(jì)算

根據(jù)文獻(xiàn)[17]思路，為避免可靠度計(jì)算對(duì)模型方差值的依賴(lài)性，采用蒙特卡羅方法進(jìn)行可靠度計(jì)算。仿真生成足夠數(shù)量的樣本退化數(shù)據(jù)，在任一時(shí)刻，定義該時(shí)刻未失效樣本數(shù)V與總樣本數(shù)N的比值為該時(shí)刻可靠度值。V是指t時(shí)刻所有邊緣退化過(guò)程的退化量都小于閾值的樣本數(shù)量。

步驟1生成N個(gè)二元性能退化過(guò)程仿真樣本。

③計(jì)算Xijk=Xi(j-1)k+ΔXijk。

④重復(fù)步驟①～步驟③N次，得到N個(gè)二元性能退化過(guò)程。

步驟2計(jì)算tj時(shí)刻可靠度值。

①比較第i個(gè)樣本退化量Xij1與D1，Xij2與D2大小。若Xij1

②對(duì)N個(gè)樣本逐次進(jìn)行步驟①，得到V最終的累加值。

③計(jì)算tj時(shí)刻可靠度值為V/N。

④重復(fù)步驟①～步驟③M次，得到所有測(cè)量時(shí)刻的可靠度值。

步驟3重復(fù)步驟1和步驟2L次，得到L個(gè)性能退化過(guò)程的可靠度曲線(xiàn)。

步驟4由L個(gè)可靠度曲線(xiàn)，進(jìn)一步可求取可靠度曲線(xiàn)的100(1-α)%置信區(qū)間。

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

電子部件通常隨導(dǎo)彈、魚(yú)水雷等裝備長(zhǎng)期貯存，在長(zhǎng)貯期間一個(gè)顯著特點(diǎn)就是產(chǎn)品個(gè)體經(jīng)歷了不同條件的貯存環(huán)境、不等次數(shù)的檢測(cè)監(jiān)測(cè)、不同時(shí)長(zhǎng)的戰(zhàn)備值班等，這些差異使得同批次電子部件在性能退化和多元相關(guān)關(guān)系等方面都呈現(xiàn)出個(gè)體差異性。此外，盡管裝備型號(hào)較多、結(jié)構(gòu)組成差異較大，但是較多的電子器件在不同電子部件的共用程度也較高，本文基于數(shù)據(jù)擬合的可靠性評(píng)估方法有較廣的適用性，因而下文將以某型裝備滾控電路板為例驗(yàn)證所提方法。

滾控電路板屬于板級(jí)電子部件，具有多個(gè)關(guān)鍵性能參數(shù)，是較為典型的多元性能退化產(chǎn)品。文獻(xiàn)[18]通過(guò)對(duì)滾控電路板失效模式分析指出負(fù)脈寬與正脈寬之差ΔT和周期T是其2個(gè)關(guān)鍵性能參數(shù)，并且由預(yù)試驗(yàn)分析結(jié)果得到加速應(yīng)力下電路板失效機(jī)理不改變的最高應(yīng)力為130℃。因而文獻(xiàn)[18]選擇9個(gè)樣本進(jìn)行恒定加速退化試驗(yàn)，確定3個(gè)加速應(yīng)力為S1=90℃，S2=103℃，S3=120℃，每個(gè)加速應(yīng)力下3個(gè)樣本。限于篇幅，詳細(xì)數(shù)據(jù)見(jiàn)文獻(xiàn)[18]。

6.1 邊緣分布參數(shù)估計(jì)

采用OpenBUGS軟件實(shí)現(xiàn)貝葉斯參數(shù)估計(jì)，所得結(jié)果如表2所示。利用邊緣退化過(guò)程隨機(jī)參數(shù)的后驗(yàn)分布樣本值，繪制箱線(xiàn)圖，如圖2所示，2個(gè)邊緣退化過(guò)程參數(shù)都具有很明顯的個(gè)體差異性，表明采用RDV模型是合理的。

表2 邊緣分布參數(shù)估計(jì)值Table 2 Parameter estimations of marginal distribution

圖2 X1和X2的隨機(jī)參數(shù)箱線(xiàn)圖Fig.2 Boxplots of random parameters of X1 and X2

6.2 隨機(jī)相關(guān)模型參數(shù)估計(jì)與模型選擇

圖3 3個(gè)加速應(yīng)力下的邊緣退化增量累積分布函數(shù)取值的散點(diǎn)圖Fig.3 Scatter plots of CDFs of degradation increments under three accelerated stresses

圖4 邊緣退化增量累積分布函數(shù)取值的散點(diǎn)圖Fig.4 Scatter plot of CDFs of degradation increments

進(jìn)一步地，分別采用Gaussian模型A，F(xiàn)rank模型A，Clayton模型A和Gumbel模型A進(jìn)行相關(guān)性建模和參數(shù)估計(jì)，參數(shù)估計(jì)值如表3所示，結(jié)果表明Frank Copula函數(shù)的DIC值最小。采用式(12)估計(jì)得到Kendallτ的非參數(shù)估計(jì)值為0.272 2，而表3同時(shí)給出了各模型τ值，可見(jiàn)Frank Copula函數(shù)τ值最為接近τ的非參數(shù)估計(jì)值。因而，不論是DIC值還是τ的非參數(shù)估計(jì)值， Frank Copula函數(shù)都是相對(duì)最優(yōu)的。

確定采用Frank Copula函數(shù)后，進(jìn)而考慮隨機(jī)相關(guān)模型。采用模型A～模型F進(jìn)行建模、參數(shù)估計(jì)，估計(jì)值如表4所示。由DIC值可以看出，模型C的DIC值最小，但是6個(gè)模型的DIC值相差都在5以?xún)?nèi)，不宜做出判斷。采用箱線(xiàn)圖進(jìn)一步判斷。圖5(a)表明，無(wú)論是模型C、模型D、模型E或者模型F，Copula函數(shù)參數(shù)都具有明顯的個(gè)體差異性；而圖5(b)～圖5(d)表明模型B、模型D、模型E和模型F，應(yīng)力變化對(duì)模型參數(shù)值影響很小。綜上，模型C是6個(gè)模型中相對(duì)最優(yōu)的模型。

表3 Copula函數(shù)參數(shù)估計(jì)值Table 3 Parameter estimations of Copula function

表4 隨機(jī)相關(guān)模型參數(shù)估計(jì)值Table 4 Parameter estimations of random correlation models

6.3 可靠性評(píng)估

額定應(yīng)力為25℃。將表2中RDV模型參數(shù)代入式(3)和式(4)，并由式(2)得到邊緣退化過(guò)程的可靠度函數(shù)R1(t)=1-FWrdv1(t)，R2(t)=1-FWrdv2(t)。將R1(t)、R2(t)及表4模型C中參數(shù)估計(jì)值代入式(8)，可得可靠度函數(shù)的解析結(jié)果Ra(t)為

圖5 隨機(jī)相關(guān)性模型參數(shù)的箱線(xiàn)圖Fig.5 Boxplots of parameters of random correlation model

(14)

式中：

按照第5節(jié)基于蒙特卡羅的算法計(jì)算可靠度值：采用模型C的參數(shù)估計(jì)值，并在區(qū)間[0,2×107]上設(shè)置測(cè)量時(shí)刻數(shù)M=4 000，單次仿真樣本數(shù)量N=200，循環(huán)次數(shù)L=200，得到仿真結(jié)果如圖6所示。圖6(a)中Rs(t)是仿真所得L條可靠度曲線(xiàn)，R(t)則是所有仿真結(jié)果Rs(t)各個(gè)時(shí)刻的可靠度均值曲線(xiàn)，Ru(t)和Rl(t)分別是仿真結(jié)果Rs(t)各個(gè)時(shí)刻的可靠度值的95%置信區(qū)間的上下限。Ra(t)是解析結(jié)果的可靠度曲線(xiàn)。盡管各個(gè)測(cè)量時(shí)刻的置信區(qū)間上下限值連線(xiàn)并不嚴(yán)格降低，但整體呈明顯下降趨勢(shì)。

圖6(a)表明，雖然解析結(jié)果在仿真結(jié)果置信區(qū)間范圍內(nèi)，但相比于仿真結(jié)果均值曲線(xiàn)，明顯低估可靠度值。對(duì)于板級(jí)電路板等高可靠性長(zhǎng)壽命產(chǎn)品的基礎(chǔ)組成部件，對(duì)其可靠度往往要求保持在0.99，甚至更高，因而計(jì)算得到了R(t)和Ra(t)可靠度值降至0.99時(shí)的壽命值分別約為194 214 h(約22年)，69 972 h(約8年),如圖6(b)所示。文獻(xiàn)[18]以測(cè)試值通常都應(yīng)服從正態(tài)分布為依據(jù)，采用了二元聯(lián)合正態(tài)分布進(jìn)行可靠性評(píng)估。可見(jiàn)文獻(xiàn)[18]對(duì)加速退化數(shù)據(jù)模型分析過(guò)于主觀(guān)，其計(jì)算得到該產(chǎn)品可靠度為0.95時(shí)的可靠壽命約為22.69年，相比于本文結(jié)果，相對(duì)保守。

圖6 可靠度曲線(xiàn)Fig.6 Reliability curves

7 結(jié) 論

本文研究了加速應(yīng)力下電子部件二元隨機(jī)相關(guān)性能退化可靠性評(píng)估，結(jié)論如下：

1) 采用考慮隨機(jī)影響的Wiener過(guò)程模型進(jìn)行邊緣分布建模，能描述性能退化過(guò)程的個(gè)體差異性，提高性能退化建模精度。Wiener過(guò)程模型隨機(jī)參數(shù)的箱線(xiàn)圖，能較好地反映性能退化過(guò)程的個(gè)體差異性，是判斷產(chǎn)品是否具有個(gè)體差異性的有效方法。

2) 產(chǎn)品個(gè)體差異性不僅僅體現(xiàn)在性能退化過(guò)程中，而且體現(xiàn)在產(chǎn)品間相關(guān)關(guān)系的差異。隨機(jī)相關(guān)模型能較好地描述相關(guān)關(guān)系的差異性，是對(duì)考慮隨機(jī)影響的隨機(jī)過(guò)程可靠性建模地進(jìn)一步完善。

3) 采用散點(diǎn)圖、DIC值和Kendallτ的非參數(shù)估計(jì)值能較好地實(shí)現(xiàn)基于Copula函數(shù)的隨機(jī)相關(guān)模型選擇。

4) 基于蒙特卡羅的可靠度計(jì)算方法，能避免可靠度計(jì)算對(duì)模型方差值的依賴(lài)，具有更好的適用性和準(zhǔn)確性。