基于Copula函數(shù)耦合性建模的二元加速退化數(shù)據(jù)統(tǒng)計(jì)分析方法

周 源，呂衛(wèi)民，王少蕾，孫 媛

(1.海軍航空大學(xué)， 山東 煙臺(tái) 264001； 2.海軍工程大學(xué)， 武漢 430000)

加速退化試驗(yàn)是一種高效評估產(chǎn)品可靠性的可行手段，已經(jīng)普遍應(yīng)用于產(chǎn)品研制、定壽、交付驗(yàn)收、延壽等多個(gè)階段[1-3]。傳統(tǒng)的加速退化試驗(yàn)只是測量產(chǎn)品某一個(gè)性能參數(shù)的退化數(shù)據(jù)，據(jù)此對產(chǎn)品可靠性進(jìn)行建模與評估。然而一些產(chǎn)品存在多種退化失效過程，并且各退化失效過程之間可能存在耦合性，產(chǎn)品失效是各退化失效過程之間競爭的結(jié)果[4-6]，因此，基于單性能參數(shù)加速退化數(shù)據(jù)的可靠性建模與評估方法已經(jīng)不適用于多參數(shù)退化產(chǎn)品。

為了利用加速退化試驗(yàn)準(zhǔn)確、高效地評估出多參數(shù)退化產(chǎn)品的可靠性，基于多元加速退化數(shù)據(jù)的統(tǒng)計(jì)分析方法逐漸受到關(guān)注[7-11]。文獻(xiàn)[7]采用多元正態(tài)分布函數(shù)擬合多元性能退化數(shù)據(jù)，在假定均值與加速應(yīng)力相關(guān)而協(xié)方差與加速應(yīng)力無關(guān)的基礎(chǔ)上建立了加速退化模型，然而，由于多元正態(tài)分布函數(shù)的適用性有限，此類方法的應(yīng)用范圍較小。文獻(xiàn)[8]采用Gamma過程與Copula函數(shù)對多元加速退化數(shù)據(jù)統(tǒng)計(jì)分析，基于多元B-S分布建立了產(chǎn)品在額定應(yīng)力下的可靠度函數(shù)。文獻(xiàn)[9]利用基于退化量分布的方法對各性能退化過程建模，假定各相同測量時(shí)刻的退化量之間具有相關(guān)性并且亦采用Copula函數(shù)進(jìn)行相關(guān)性建模。

繼電器廣泛用于國防工業(yè)的多種設(shè)備，主要功能是可靠地完成電信號傳遞、電路隔離、繼電保護(hù)等，繼電器的可靠性對整個(gè)設(shè)備能否正常工作非常重要[12-13]。對于很多用于軍工、航天的繼電器產(chǎn)品，會(huì)有一段比較長的貯存期，繼電器在長期貯存過程中的某些性能指標(biāo)不可避免會(huì)產(chǎn)生退化[13-15]。因此，如何對繼電器的性能參數(shù)退化量進(jìn)行合理統(tǒng)計(jì)分析，從而掌握繼電器的質(zhì)量變化規(guī)律是非常必要的。目前的性能退化數(shù)據(jù)統(tǒng)計(jì)分析方法大都是針對一種性能退化參數(shù)，然而工程經(jīng)驗(yàn)表明繼電器具有明顯退化趨勢的退化參數(shù)包括接觸電阻值、釋放電壓。此外，由于繼電器內(nèi)部的結(jié)構(gòu)相關(guān)性，接觸電阻值、釋放電壓之間具有一定的耦合性，即任一種性能參數(shù)的退化過程會(huì)受到其他性能參數(shù)退化過程的影響。本文提出一種基于Copula函數(shù)耦合性建模的二元加速退化數(shù)據(jù)統(tǒng)計(jì)分析方法，結(jié)合加速退化試驗(yàn)評估某型繼電器產(chǎn)品的可靠性。

1 二元加速退化數(shù)據(jù)建模

某型繼電器存在兩種性能退化過程，令Xp,k(t)表示產(chǎn)品在溫度應(yīng)力水平Tk下的第p個(gè)性能退化過程，其中p=1,2，任一性能退化過程首次達(dá)到失效閾值Dp時(shí)產(chǎn)品發(fā)生退化失效。根據(jù)前期的研究結(jié)論，Xp,k(t)服從Wiener隨機(jī)過程，并且X1,k(t)與X2,k(t)之間具有耦合性。Tk代表第k個(gè)加速溫度應(yīng)力水平，ti代表第i次測量時(shí)間，xpjk(ti)表示Tk下第j個(gè)樣品的第p個(gè)性能退化過程Xp,k(t)在時(shí)間ti的退化測量值，其中k=1,2,3；j=1,2,…,Nk；p=1,2；i=1,2,…,Mk。Tk下所有樣品的性能退化測量值利用如下的數(shù)據(jù)結(jié)構(gòu)描述：

(1)

(2)

(3)

(4)

式(4)中，Φ(·)為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)。根據(jù)文獻(xiàn)[2]中的研究結(jié)論，Wiener過程的漂移參數(shù)值與擴(kuò)散參數(shù)值都與環(huán)境應(yīng)力水平相關(guān)，并且可利用加速模型定量表示它們的關(guān)系，如果環(huán)境應(yīng)力為絕對溫度Tk，可采用阿倫尼斯加速模型表示為

μpk=exp(ap-bp/Tk)

(5)

(6)

式(5)、式(6)中，ap,bp,cp為第p個(gè)性能參數(shù)對應(yīng)的加速模型待估參數(shù)。

設(shè)G(x1)，G(x2)分別為變量X1，X2的累積分布函數(shù)，根據(jù)Sklar’s理論，如果H(x1,x2)為G(x1)，G(x2)的聯(lián)合分布函數(shù)，則存在一個(gè)滿足下式的Copula函數(shù)C(·)，

H(x1,x2)=C(G(x1),G(x2);θ)

(7)

式(7)中，θ為Copula參數(shù)。目前應(yīng)用較廣的Copula函數(shù)為Gaussian，F(xiàn)rank，Gumbel，Clayton等，這4種Copula函數(shù)具有不同的耦合效果[16-18]，表1給出了這4種Copula的分布函數(shù)及密度函數(shù)。

H(x1,x2)的概率密度函數(shù)為

c(G(x1),G(x2))g(x1)g(x2)

(8)

將式(4)，式(5)，式(6)代入式(7)，得

H(Δx1jk(ti),Δx2jk(ti))=

C(Φ(U1jk(ti)),Φ(U2jk(ti));θk)

(9)

將式(3)，式(9)代入式(7)，得

h(Δx1jk(ti),Δx2jk(ti))=

c(Φ(U1jk(ti)),Φ(U2jk(ti));θk)g(Δx1jk(ti))g(Δx2jk(ti))

(10)

2 參數(shù)估計(jì)

參數(shù)估計(jì)的任務(wù)是提出有效的方法估計(jì)出二元加速退化模型的各待估參數(shù)值，包括以下10個(gè)參數(shù)：ap，bp，cp，Λp，e，f，其中p=1,2。根據(jù)式(10)建立如下極大似然方程：

g(Δx1jk(ti))g(Δx2jk(ti))}

(11)

然而由于待估參數(shù)過多，難以通過傳統(tǒng)的極大似然估計(jì)法直接極大化式(11)，無法一體化獲取參數(shù)估計(jì)值。因此考慮采用Bayesian MCMC的思路[19]，設(shè)計(jì)如下的參數(shù)估計(jì)算法一體化估計(jì)出各參數(shù)值，主要步驟為：

1) 令logLijk=log{c(Φ(U1jk(ti)),Φ(U2jk(ti));θk)g(Δx1jk(ti))g(Δx2jk(ti))}；

2) 構(gòu)建服從Poisson分布的偽隨機(jī)變量Yijk，設(shè)Yijk的所有觀測值為0并且Poisson分布的均值參數(shù)為-logLijk+R，如Yijk～Poi(0;-logLijk+R)，R為一個(gè)滿足-logLijk+R>0的實(shí)數(shù)；

3) 由Yijk的概率密度函數(shù)fP(0;-logLijk+R)建立如下極大似然方程，

(12)

4) 根據(jù)以下關(guān)系式將Poisson分布參數(shù)估計(jì)問題等效為二元加速退化模型參數(shù)估計(jì)問題

(13)

5) 在OPENBUGS軟件中建立基于Poisson分布Yijk～Poi(0;-logLijk+R)的貝葉斯模型，將ap，bp，cp，Λp，e，f分別表示為服從無信息先驗(yàn)分布的隨機(jī)參數(shù)，并設(shè)定各待估參數(shù)的初值；

6) 利用基于Gibbs抽樣的馬爾科夫鏈蒙特卡洛方法經(jīng)過足夠多步的迭代，獲得ap，bp，cp，Λp，e的后驗(yàn)均值；

7) 判斷各隨機(jī)參數(shù)的后驗(yàn)分布是否收斂，如果收斂則將后驗(yàn)均值作為參數(shù)估計(jì)值，否則調(diào)整各隨機(jī)參數(shù)的先驗(yàn)分布、待估參數(shù)初值、迭代步數(shù)，直到各隨機(jī)參數(shù)的后驗(yàn)分布都達(dá)到收斂。

3 可靠性預(yù)計(jì)

(14)

R0(t)=P(X1，0(t)

(15)

式(15)中，

(16)

如果X1,k(t)與X2,k(t)之間沒有耦合性，則繼電器的可靠度函數(shù)為

(17)

4 案例應(yīng)用

為了掌握某型繼電器在長期貯存過程中的可靠性變化規(guī)律，首先結(jié)合以往的工程經(jīng)驗(yàn)分析了此型繼電器產(chǎn)品的失效模式與失效機(jī)理，得出：接觸電阻增大、釋放電壓降低是產(chǎn)品在貯存過程中最主要的兩種退化失效模式；接觸電阻增大與釋放電壓降低兩種失效過程存在某種程度的耦合性；溫度是產(chǎn)品接觸電阻增大、釋放電壓降低的敏感環(huán)境應(yīng)力。然后，隨機(jī)選取了24個(gè)樣品開展了高溫加速退化試驗(yàn)。具體試驗(yàn)信息為：在正常貯存溫度25 ℃(298.16 K)下分別測量每個(gè)樣品的接觸電阻初值x0、釋放電壓初值y0；24個(gè)樣品被平均分配到3個(gè)老化試驗(yàn)箱內(nèi)，箱內(nèi)的溫度分別設(shè)為60 ℃(333.16 K)，90 ℃(363.16 K)，120 ℃(393.16 K)；試驗(yàn)過程中每隔0.4×103h測量一次所有樣品的接觸電阻值x與釋放電壓值y，試驗(yàn)截止時(shí)間為2.8×103h。

圖1、圖2、圖3展示了樣品在3個(gè)溫度應(yīng)力水平下的接觸電阻值與釋放電壓值的變化曲線。當(dāng)接觸電阻測量值相對于初值的變化量首次達(dá)到20 mΩ時(shí)或當(dāng)釋放電壓測量值相對于初值的變化量首次達(dá)到2 V時(shí)，產(chǎn)品會(huì)發(fā)生退化失效。以下利用本文所提方法對二元加速退化數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，預(yù)測此型繼電器產(chǎn)品在正常貯存溫度下的貯存可靠度變化規(guī)律。

圖1 繼電器在333.16 K下的接觸電阻及釋放電壓退化數(shù)據(jù)

圖2 繼電器在363.16 K下的接觸電阻及釋放電壓退化數(shù)據(jù)

圖3 繼電器在393.16 K下的接觸電阻及釋放電壓退化數(shù)據(jù)

首先，利用第1節(jié)中所提方法對此型繼電器產(chǎn)品建立二元耦合加速退化模型，如式(9)所示，其中Copula函數(shù)類型分別選用Clayton Copula、Frank Copula。然后，利用第2節(jié)中所提方法進(jìn)行參數(shù)估計(jì)。當(dāng)選用Clayton Copula時(shí)，估計(jì)出各參數(shù)值如表2所示，當(dāng)選用Frank Copula時(shí)，估計(jì)出各參數(shù)值如表3所示，其中Mean為隨機(jī)參數(shù)的后驗(yàn)估計(jì)均值，SD為后驗(yàn)估計(jì)標(biāo)準(zhǔn)差，MC_error為MCMC參數(shù)估計(jì)的誤差，2.5%表示后驗(yàn)估計(jì)值的2.5%分位數(shù),97.5%表示后驗(yàn)估計(jì)值的97.5%分位數(shù)。比較表2與表3中的參數(shù)估計(jì)結(jié)果，各參數(shù)的后驗(yàn)均值雖然相差不大但均有一定變化，說明選用不同的Copula函數(shù)會(huì)影響最終的參數(shù)估計(jì)結(jié)果。通過比較AIC(Akaike Information Criterion)值，選用Clayton Copula的二元加速退化模型能夠與試驗(yàn)數(shù)據(jù)擬合得更優(yōu)。

表2 選用Clayton Copula時(shí)的參數(shù)估計(jì)

表3 選用Frank Copula時(shí)的參數(shù)估計(jì)

上式中，

表4 接觸電阻加速退化模型的參數(shù)估計(jì)

表5 釋放電壓加速退化模型的參數(shù)估計(jì)

圖4 可靠度評估結(jié)果

5 結(jié)論

1) 針對某型繼電器產(chǎn)品，提出了一種可行的二元加速退化數(shù)據(jù)建模方法，包括利用Wiener過程建立性能參數(shù)的退化模型；結(jié)合阿倫尼斯方程建立模型參數(shù)的加速退化模型；采用Copula函數(shù)建立二元加速退化過程之間的耦合性模型。

2) 提出了一種基于Bayesian MCMC的參數(shù)估計(jì)，有效解決了因?yàn)槎铀偻嘶Ｐ椭写绤?shù)過多，傳統(tǒng)估計(jì)方法不適用的難題。

3) Copula函數(shù)能夠描述兩種退化過程之間的耦合性，但也給可靠性建模與參數(shù)估計(jì)增加了難度。

4) 通過所提方法評估了某型繼電器的可靠性，與工程實(shí)際中的認(rèn)識相符，但評估結(jié)果還需產(chǎn)品的真實(shí)貯存失效數(shù)據(jù)做進(jìn)一步驗(yàn)證。

[1] LI J X,WANG Z H,LIU X,et al.A Wiener process model for accelerated degradation analysis considering measurement errors[J].Microelectronics Reliability,2016,65(10):8-15.

[2] 滕飛,王浩偉,陳瑜,等.加速度計(jì)加速退化數(shù)據(jù)統(tǒng)計(jì)分析方法[J].中國慣性技術(shù)學(xué)報(bào),2017,25(2):275-280.

[3] 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.

[4] 王浩偉,奚文駿,馮玉光.基于退化失效與突發(fā)失效競爭的導(dǎo)彈剩余壽命預(yù)測[J].航空學(xué)報(bào),2016,37(4):1240-1248.

[5] 王浩偉,滕克難,李軍亮.隨機(jī)環(huán)境應(yīng)力沖擊下基于多參數(shù)相關(guān)退化的導(dǎo)彈部件壽命預(yù)測[J].航空學(xué)報(bào),2016,37(11):3404-3412.

[6] WANG Y P,PHAM H.Modeling the dependent competing risks with multiple degradation processes and random shock using time-varying Copulas[J].IEEE Transactions on Reliability,2012,61(1):13-22.

[7] 潘駿,王小云,陳文華,等.基于多元性能參數(shù)的加速退化試驗(yàn)方案優(yōu)化設(shè)計(jì)研究[J].機(jī)械工程學(xué)報(bào),2012,48(2):30-35.

[8] PAN Z Q,BALAKRISHNAN N.Reliability modeling of degradation of products with multiple performance characteristics based on gamma process[J].Reliability Engineering and System Safety,2011,96(8):949-957.

[9] 王孟.基于Copula函數(shù)的多元加速退化試驗(yàn)方法研究[D].杭州：浙江理工大學(xué),2013.

[10] TSAI T R,SUNG W Y,LIO Y L,et al.Optimal two-variable accelerated degradation test plan for Gamma degradation processes[J].IEEE Transactions on Reliability,2016,65(1):459-468.

[11] ZHANG X P,SHANG J Z,CHEN X,et al.Statistical inference of accelerated life testing with dependent competing failures based on Copula theory[J].IEEE Transactions on Reliability,2014,63(3):764-780.

[12] 李文華，馬思寧，沈培根，等.振動(dòng)條件下鐵路繼電器壽命預(yù)測研究[J].電氣工程學(xué)報(bào),2017,12(7):8-15.

[13] 葉雪榮，林義剛，黃曉毅，等.航天繼電器貯存過程吸合時(shí)間退化機(jī)理研究[J].電工技術(shù)學(xué)報(bào),2017,32(11):173-179.

[14] 孟彥辰.基于性能參數(shù)退化的航天繼電器可靠性建模與評估研究[D].哈爾濱：哈爾濱工業(yè)大學(xué),2015.

[15] 王浩偉,滕克難.基于加速退化數(shù)據(jù)的可靠性評估技術(shù)綜述[J].系統(tǒng)工程與電子技術(shù),2017,39(12):2877-2885.

[16] FAN M F,ZENG Z G,ZIO E,et al.Modeling dependent competing failure processes with degradation-shock dependence[J].Reliability Engineering and System Safety,2017,165(7):422-430.

[17] 潘正強(qiáng).加速應(yīng)力下二元退化可靠性建模及其試驗(yàn)設(shè)計(jì)方法[D].長沙:國防科學(xué)技術(shù)大學(xué),2011.

[18] PENG W W,LI Y F,YANG Y J,et al.Bivariate analysis of incomplete degradation observations based on inverse Gaussian processes and Copulas[J].IEEE Transactions on Reliability,2016,65(2):624-639.

[19] NTZOUFRAS I.Bayesian Modeling Using WinBUGS[M].Hoboken,New Jersey,John Wiley & Sons,2009.

[20] PAN Z Q,BALAKRISHNAN N,SUN Q.Bivariate constant-stress accelerated degradation model and inference[J].Communications in Statistics-Simulation and Computation,2011,40(2):247-257.