李 玥, 冶繼民

(西安電子科技大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 西安 710126)

在某些實(shí)際問題中, 為了解系統(tǒng)的失效原因和隨時(shí)間變化系統(tǒng)的退化情況, 需要在試驗(yàn)的不同階段將部分未失效系統(tǒng)移離試驗(yàn), 即逐步截尾試驗(yàn). 逐步截尾方案通常分為逐步Ⅰ型截尾方案(progressively Type-Ⅰ censored scheme, Type-Ⅰ PCS)和逐步Ⅱ型截尾方案(progressively Type-Ⅱ censored scheme, Type-Ⅱ PCS). 逐步Ⅰ型截尾是指在預(yù)先給定的間隔觀察時(shí)刻點(diǎn), 將預(yù)先給定數(shù)量的未失效系統(tǒng)移離試驗(yàn), 直到試驗(yàn)終止時(shí)刻結(jié)束試驗(yàn). 逐步Ⅱ型截尾是指每當(dāng)觀測(cè)到系統(tǒng)失效時(shí), 將預(yù)先給定數(shù)量的未失效系統(tǒng)移離試驗(yàn), 直到有預(yù)先給定數(shù)量的系統(tǒng)失效時(shí)結(jié)束試驗(yàn). 逐步截尾方案比傳統(tǒng)的Ⅰ型和Ⅱ型截尾方案更靈活, 且Type-Ⅱ PCS與Type-Ⅰ PCS相比, 可以保證至少觀測(cè)到給定數(shù)量的失效系統(tǒng), 從而提高試驗(yàn)的準(zhǔn)確性[1-2].

系統(tǒng)壽命數(shù)據(jù)包括失效時(shí)間和失效原因兩方面. 在工程實(shí)踐中, 引起系統(tǒng)失效的原因多種多樣, 但由于測(cè)試環(huán)境、 測(cè)試時(shí)間等因素的限制, 有時(shí)可能無(wú)法識(shí)別導(dǎo)致系統(tǒng)失效的確切部件, 而只能將其歸咎于一些可能導(dǎo)致系統(tǒng)失效的部件集合. 此時(shí), 系統(tǒng)的失效原因被屏蔽, 這種數(shù)據(jù)稱為屏蔽數(shù)據(jù), 相應(yīng)的系統(tǒng)稱為屏蔽數(shù)據(jù)系統(tǒng). Usher等[3]采用極大似然方法研究了串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的可靠性; Sarhan[4-5]研究了當(dāng)部件壽命服從Weibull分布和Pareto分布時(shí), 屏蔽數(shù)據(jù)系統(tǒng)模型的統(tǒng)計(jì)推斷; 蔡靜等[6]研究了Burr Ⅻ分布串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的可靠性; Sarhan等[7]討論了并聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的部件可靠性評(píng)估問題.

上述研究均假設(shè)系統(tǒng)部件壽命變量相互獨(dú)立, 但在工程實(shí)踐中, 系統(tǒng)部件之間通常不獨(dú)立, 部件失效存在一定的相關(guān)關(guān)系. 例如, 一架雙引擎飛機(jī), 當(dāng)其中一個(gè)引擎出現(xiàn)故障時(shí), 另一個(gè)引擎的工作狀態(tài)會(huì)發(fā)生變化, 從而維持飛機(jī)能正常飛行. 因此, 研究部件之間的相依關(guān)系及其對(duì)系統(tǒng)可靠性的影響有一定意義. Feizjavadian等[8]和Bai等[9]利用二元Marshall-Olkin Weibull分布描述失效原因之間的相依結(jié)構(gòu), 對(duì)相依串聯(lián)系統(tǒng)模型進(jìn)行了統(tǒng)計(jì)分析； Cai等[10]基于二元Marshall-Olkin Weibull分布對(duì)相依串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)進(jìn)行了統(tǒng)計(jì)分析； 付倩嬈等[11]討論了二元Marshall-Olkin指數(shù)分布下含有屏蔽數(shù)據(jù)的相依串聯(lián)系統(tǒng)的可靠性分析問題. 上述研究均采用多元分布刻畫相依結(jié)構(gòu), 適用于系統(tǒng)部件壽命變量具有相同分布的結(jié)構(gòu). 但在復(fù)雜結(jié)構(gòu)的情形下, 多元分布對(duì)變量之間的不對(duì)稱性和尾部相依性的計(jì)算較困難. Copula函數(shù)因其靈活性可以很好地解決上述問題, 本文選擇Copula函數(shù)進(jìn)行相依結(jié)構(gòu)的建模. Hsu等[12]基于Clayton Copula函數(shù)研究了在Ⅰ型截尾試驗(yàn)下多部件相依串聯(lián)系統(tǒng)的可靠性評(píng)估問題； 蔡靜等[13]基于Clayton Copula函數(shù)討論了相依并-串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的可靠性分析問題.

Burr Ⅻ分布[14]是一種重要的連續(xù)壽命分布. 由于其失效率函數(shù)的非單調(diào)形式, 使其適合代表許多產(chǎn)品的壽命, 具有很好的靈活性, 廣泛應(yīng)用于可靠性工程、 通信工程、 質(zhì)量控制、 航空航天等領(lǐng)域. Burr Ⅻ分布作為一種失效模型, 關(guān)于其應(yīng)用和統(tǒng)計(jì)性質(zhì)目前已有許多研究結(jié)果. 但基于Burr Ⅻ分布對(duì)相依屏蔽數(shù)據(jù)系統(tǒng)的研究報(bào)道較少, 因此, 本文針對(duì)部件壽命變量服從Burr Ⅻ分布的情形, 通過(guò)引入Copula函數(shù)建立部件壽命變量之間的相依結(jié)構(gòu), 并討論在逐步Ⅱ型截尾試驗(yàn)下相依串聯(lián)屏蔽數(shù)據(jù)系統(tǒng)的可靠性. 首先, 推導(dǎo)出模型參數(shù)及系統(tǒng)可靠度函數(shù)的極大似然估計(jì)(maximum likelihood estimations, MLEs)； 然后, 基于漸近正態(tài)性理論和bootstrap抽樣算法分別構(gòu)造出參數(shù)的漸近置信區(qū)間(asymptotic confidence intervals, ACIs)及偏差校正的百分位bootstrap置信區(qū)間(bias-corrected percentile bootstrap confidence intervals, Boot-BCP CIs)； 最后, 用仿真模擬和真實(shí)數(shù)據(jù)分析驗(yàn)證所提出統(tǒng)計(jì)方法的可行性和有效性.

1 Copula理論

1.1 Copula函數(shù)

Copula函數(shù)是描述多元隨機(jī)變量間相依結(jié)構(gòu)的一種函數(shù), 它可以將多元隨機(jī)變量的聯(lián)合分布函數(shù)與其邊際分布函數(shù)聯(lián)系起來(lái). 下面簡(jiǎn)要介紹二元變量情形下Copula函數(shù)的基本理論[15-16].

失效Copula函數(shù)是一個(gè)多維聯(lián)合分布函數(shù), 其定義域?yàn)閇0,1]×…×[0,1]. 根據(jù)Sklar定理[15], 它可以連接多個(gè)隨機(jī)變量的邊際分布構(gòu)造聯(lián)合分布. 對(duì)于二元變量的情形, 其可表示為

H(z1,z2)=C(F1(z1),F2(z2)),

(1)

二元生存Copula函數(shù)和二元失效Copula函數(shù)的關(guān)系為

1.2 相依性度量

通常的相依性度量如Spearman’sρ和Kendall’sτ都可采用Copula函數(shù)表示[15]為

1.3 Archimedean Copula族

由于Archimedean Copula族可以準(zhǔn)確地描述非對(duì)稱隨機(jī)變量之間的尾相關(guān)性, 因此在實(shí)際生活中應(yīng)用廣泛. 二元Archimedean Copula族函數(shù)定義[15]為C(μ,ν)=φ-1(φ(μ)+φ(ν)), 其中生成函數(shù)φ(·)是一個(gè)嚴(yán)格單調(diào)遞減的凸函數(shù), 滿足φ(μ)+φ(ν)≤φ(0),φ(1)=0,μ,ν∈[0,1],φ-1(·)是φ(·)的逆.

Clayton Copula隸屬于Archimedean Copula族, 其Copula函數(shù)為

CC(μ,ν)=[μ-θ+ν-θ-1]-1/θ,θ∈[-1,∞){0}.

Clayton Copula適用于描述下尾相關(guān)、 上尾漸近獨(dú)立隨機(jī)變量之間的相關(guān)性, 由Kendall’sτ及相應(yīng)Copula函數(shù)的關(guān)系可得τ=θ/(θ+2), 當(dāng)θ→0時(shí), 隨機(jī)變量趨于相互獨(dú)立.

2 模型描述和基本假設(shè)

2.1 模型描述

由于某些系統(tǒng)的失效原因并不確切, 因此用Ki(i=1,2,…,m)表示系統(tǒng)i可能失效原因的集合,ki是Ki的觀測(cè)值.ki={1}表示系統(tǒng)i失效是由部件1導(dǎo)致的；ki={2}表示系統(tǒng)i失效是由部件2導(dǎo)致的；ki={1,2}表示系統(tǒng)i的失效原因被屏蔽.因此, 觀測(cè)到的系統(tǒng)壽命數(shù)據(jù)為(t1,k1,r1),(t2,k2,r2),…,(tm,km,rm).

2.2 基本假設(shè)

1) 設(shè)Z1,Z2分別表示部件1和2的失效時(shí)間, 各部件的失效時(shí)間Zv(v=1,2)服從參數(shù)為γv,α的Burr Ⅻ分布, 簡(jiǎn)記為Zv～BX(γv,α), 其中γv>0,α>0.系統(tǒng)的失效時(shí)間為各部件失效時(shí)間的最小值, 即T=min{Z1,Z2}.

部件v失效時(shí)間的累積分布函數(shù)(cumulative distribution function, CDF)、 概率密度函數(shù)(probability density function, PDF)、 生存函數(shù)(survival function, SF)(也稱可靠度函數(shù))分別為

2) 假設(shè)系統(tǒng)由兩個(gè)相依的部件串聯(lián)而成, 且其失效是由最先失效的部件所致.為考察數(shù)據(jù)間的下尾相依性, 本文選擇Clayton Copula刻畫各部件之間的相依關(guān)系, 則系統(tǒng)的SF為

S(t)=SZ1,Z2(z1,z2)|(t,t)=P(Z1>t,Z2>t)=[(1+tα)γ1θ+(1+tα)γ2θ-1]-1/θ.

3) 等概率假設(shè)： 屏蔽的發(fā)生與失效原因無(wú)關(guān), 即

P(Ki=ki|Ti=ti,Vi=v′)=P(Ki=ki|Ti=ti,Vi=v), ?v′,v∈{1,2},

其中,Ti和Vi分別表示系統(tǒng)i的失效時(shí)間和確切失效原因,p?P(Ki=ki|Ti=ti,Vi=v)表示屏蔽概率.為計(jì)算簡(jiǎn)便, 本文基于等概率假設(shè)進(jìn)行討論, 也有部分研究放寬了該假設(shè)[12,17-18].

3 極大似然估計(jì)

由部件1導(dǎo)致系統(tǒng)i失效的PDF為

由部件2導(dǎo)致系統(tǒng)i失效的PDF為

則失效系統(tǒng)的似然函數(shù)可表示為

(4)

為簡(jiǎn)化似然函數(shù), 定義失效原因示性函數(shù)為

定義輔助示性函數(shù)為

因此基于假設(shè)3), 似然函數(shù)可表示為

對(duì)數(shù)似然函數(shù)為

對(duì)lnL(γ1,γ2,α,θ)關(guān)于γ1,γ2,α,θ求一階導(dǎo)數(shù)并令其等于0, 可得似然方程為

由于很難得到模型參數(shù)的顯式表達(dá)式, 因此本文考慮迭代算法, 如Newton-Raphson迭代算法或其他迭代算法求解γ1,γ2,α,θ的MLEs.則系統(tǒng)可靠度函數(shù)的估計(jì)可表示為

(11)

4 區(qū)間估計(jì)

4.1 漸近置信區(qū)間

記參數(shù)γ1,γ2,α,θ的Fisher信息矩陣為

其中,L(·,·)表示lnL(γ1,γ2,α,θ)對(duì)參數(shù)γ1,γ2,α,θ的二階偏導(dǎo)數(shù), 其中,

則參數(shù)γ1,γ2,α,θ的100(1-ε)%雙側(cè)ACIs為

其中Zε/2是標(biāo)準(zhǔn)正態(tài)分布的ε/2分位數(shù).

4.2 Bootstrap置信區(qū)間

在現(xiàn)代制造業(yè)中, 由于時(shí)間和成本的限制, 樣本量通常不大, 基于大樣本數(shù)據(jù)的統(tǒng)計(jì)方法可能并不合適, 有時(shí)甚至?xí)a(chǎn)生誤導(dǎo). 因此, 本文采用bootstrap抽樣算法構(gòu)造參數(shù)的置信區(qū)間.

5) 參數(shù)γ1,γ2,α,θ的100(1-ε)%雙側(cè)Boot-BCP CIs為

φ-1(·)是φ(·)的逆函數(shù),I(·)表示示性函數(shù).

5 仿真模擬

下面利用Monte Carlo仿真模擬評(píng)估上述統(tǒng)計(jì)方法的性能并分析模擬結(jié)果. 考慮一個(gè)具有兩個(gè)相依部件的串聯(lián)系統(tǒng), 假設(shè)部件的壽命分布服從Burr Ⅻ分布, 利用二維Clayton Copula函數(shù)刻畫其相依性. 選取參數(shù)γ1=1,γ2=1.2,α=1, 給定樣本量n和失效系統(tǒng)數(shù)m以及截尾方案{r1=r2=…=rm-1=0,rm=n-m}, 在不同的屏蔽概率(p=0.3,0.5)和不同的相依關(guān)系(θ=2,4)下, 重復(fù)試驗(yàn)5 000次, 并計(jì)算各參數(shù)的MLEs、 均方誤差(mean square errors, MSEs), 以及ACIs和Boot-BCP CIs的平均寬度(average widths, AWs)、 覆蓋率(coverage probabilities, CPs), 結(jié)果分別列于表1～表3. 當(dāng)時(shí)間t=1,2,3時(shí)系統(tǒng)的可靠度分別為0.347 6,0.213 4,0.152 8.

表1 當(dāng)p=0.3, θ=2時(shí), 不同樣本量下模型參數(shù)的點(diǎn)估計(jì)和區(qū)間估計(jì)

表2 當(dāng)p=0.5, θ=2時(shí), 不同樣本量下模型參數(shù)的點(diǎn)估計(jì)和區(qū)間估計(jì)

表3 當(dāng)p=0.3, θ=4時(shí), 不同樣本量下模型參數(shù)的點(diǎn)估計(jì)和區(qū)間估計(jì)

由表1～表3可見: 隨著樣本量的增大, 點(diǎn)估計(jì)中參數(shù)的MLEs逐漸接近真值, MSEs逐漸減?。?區(qū)間估計(jì)中, ACI-AWs和BCI-AWs逐漸減小, ACI-CPs和BCI-CPs逐漸增大, 這是因?yàn)殡S著n或m的增加, 試驗(yàn)得到的有效樣本量增多, 從而估計(jì)結(jié)果更準(zhǔn)確. 隨著屏蔽概率的增大, MSEs,AWs和CPs的結(jié)果表明, 參數(shù)點(diǎn)估計(jì)和區(qū)間估計(jì)的準(zhǔn)確性都會(huì)降低, 這是因?yàn)槠帘螖?shù)據(jù)越多, 試驗(yàn)得到的有效信息越少, 因此在分析壽命數(shù)據(jù)時(shí), 屏蔽數(shù)據(jù)的存在不可忽略. 隨著部件壽命變量之間相依性的增強(qiáng), 參數(shù)的MLEs接近真值, ACI-AWs,BCI-AWs都減小, ACI-CPs,BCI-CPs都增大, 表明在一定程度上考慮部件壽命變量之間的相依性是必要的.

為進(jìn)一步說(shuō)明屏蔽概率對(duì)系統(tǒng)可靠度估計(jì)的影響, 表4列出了當(dāng)p=0.3和p=0.5時(shí)系統(tǒng)可靠度估計(jì)值以及其與真實(shí)值的誤差.由表4可見, 隨著樣本量的增加, 系統(tǒng)可靠度估計(jì)值越來(lái)越接近真值.屏蔽概率的增大會(huì)導(dǎo)致降低系統(tǒng)可靠度估計(jì)的精度.

表4 當(dāng)p=0.3和p=0.5時(shí)系統(tǒng)可靠度的比較

當(dāng)p=0.5,θ=2時(shí), 刪除失效原因被屏蔽的數(shù)據(jù), 基于失效部件明確的數(shù)據(jù)構(gòu)造參數(shù)的極大似然估計(jì), 所得結(jié)果列于表5. 對(duì)比表2和表5可見, 參數(shù)極大似然估計(jì)的MSE增大. 表明忽略屏蔽數(shù)據(jù)將導(dǎo)致增大參數(shù)估計(jì)的誤差, 因此, 屏蔽數(shù)據(jù)的存在不可忽略.

表5 當(dāng)p=0.5, θ=2時(shí), 刪除屏蔽數(shù)據(jù)后模型參數(shù)的估計(jì)

6 實(shí)例分析

考慮用文獻(xiàn)[21]中對(duì)36個(gè)小型電器進(jìn)行壽命試驗(yàn)的實(shí)際壽命數(shù)據(jù)集驗(yàn)證本文方法的可行性. 在該數(shù)據(jù)集中, 失效原因被分為18種不同的模式, 觀察到的大多數(shù)故障都是由模式9引起的, 因此本文將失效模式9視為失效原因1, 將其他所有失效模式視為失效原因2. 給定n=36,m=30以及截尾方案{r1=6,r2=r3=…=r30=0}, 得到含有屏蔽數(shù)據(jù)的逐步Ⅱ型截尾數(shù)據(jù)列于表6.

表6 含有屏蔽數(shù)據(jù)的逐步Ⅱ型截尾數(shù)據(jù)

基于含有屏蔽數(shù)據(jù)的逐步Ⅱ型截尾數(shù)據(jù), 模型參數(shù)γ1,γ2,α,θ相應(yīng)的MLEs,ACIs和Boot-BCP CIs列于表7. 由表7可見, ACIs和Boot-BCP CIs均包含相應(yīng)的點(diǎn)估計(jì), Boot-BCP CIs在區(qū)間長(zhǎng)度上性能優(yōu)于ACIs.

表7 模型參數(shù)估計(jì)

為檢驗(yàn)本文模型的擬合優(yōu)度, 采用Kolmogorov-Smirnov(K-S)檢驗(yàn)的K-S距離和p值(表8)以及如圖1所示的Quantile-Quantile(Q-Q)圖說(shuō)明Burr Ⅻ分布的擬合情況. 由表8可見, 當(dāng)選取置信度為0.05時(shí), 原因1和原因2的K-S距離均小于相應(yīng)的臨界值, 且相應(yīng)的p值均大于0.05. Q-Q圖也表明, 擬合的Burr Ⅻ分布是一個(gè)合適的分布. 因此, 本文提出的模型可行且有效.

表8 K-S檢驗(yàn)

圖1 Q-Q圖Fig.1 Q-Q plots

表9 Copula模型與獨(dú)立模型的比較

綜上所述, 本文針對(duì)部件壽命變量服從Burr Ⅻ分布的情形, 基于Clayton Copula函數(shù), 討論了含有屏蔽數(shù)據(jù)的相依串聯(lián)系統(tǒng)的可靠性, 得到了模型參數(shù)及系統(tǒng)可靠度函數(shù)的MLEs, 并基于漸近正態(tài)性理論和bootstrap抽樣算法構(gòu)造了模型參數(shù)的ACIs和Boot-BCP CIs. 仿真模擬和真實(shí)數(shù)據(jù)分析驗(yàn)證了所提出統(tǒng)計(jì)方法的可行性與有效性. 在小樣本情形下, 建議采用偏差校正的百分位bootstrap方法進(jìn)行區(qū)間估計(jì). 在系統(tǒng)壽命試驗(yàn)的分析中屏蔽數(shù)據(jù)的存在不可忽略, 隨著屏蔽概率的增大模型參數(shù)估計(jì)的精度逐漸降低.