黃首清，劉慶海*，李芳勇，楊 勇，遇 今，張兆霖

（1.航天機(jī)電產(chǎn)品環(huán)境可靠性試驗(yàn)技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室； 2.可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室；3.北京衛(wèi)星環(huán)境工程研究所：北京 100094）

0 引言

單應(yīng)力多個(gè)應(yīng)力水平加速壽命試驗(yàn)是經(jīng)典的加速壽命試驗(yàn)之一；其可根據(jù)所設(shè)計(jì)的多個(gè)應(yīng)力水平的加速壽命試驗(yàn)失效時(shí)間數(shù)據(jù)，在失效機(jī)理不變的前提下外推其他應(yīng)力水平的壽命。單應(yīng)力多個(gè)應(yīng)力水平加速壽命試驗(yàn)涉及到應(yīng)力水平的數(shù)量和大小的選取、各應(yīng)力水平壽命分布擬合、失效機(jī)理一致性判別、加速模型的選用和擬合、加速模型與壽命分布模型的關(guān)聯(lián)，以及壽命外推和可靠壽命計(jì)算等多個(gè)環(huán)節(jié)，對(duì)相關(guān)技術(shù)和工程經(jīng)驗(yàn)要求較高，因此，簡化應(yīng)用復(fù)雜度以及開展軟件設(shè)計(jì)研究是推廣該技術(shù)的重要途徑和研究方向。

在單應(yīng)力多個(gè)應(yīng)力水平加速壽命試驗(yàn)領(lǐng)域的相關(guān)研究主要聚焦在試驗(yàn)的精細(xì)化設(shè)計(jì)和數(shù)據(jù)的高精度分析[1-2]。早期的研究集中在應(yīng)力水平和樣本量的選擇。例如：Meeker 等[3]認(rèn)為對(duì)于3 個(gè)應(yīng)力水平的情況，從低應(yīng)力水平到高應(yīng)力水平的樣本量比例應(yīng)為4:2:1，且中等應(yīng)力水平應(yīng)為低應(yīng)力水平和高應(yīng)力水平的平均值；而中國國家標(biāo)準(zhǔn)GB 2689.1—1981 [4]規(guī)定，在缺少加速因子相關(guān)歷史數(shù)據(jù)的情況下，需要設(shè)計(jì)多個(gè)應(yīng)力水平的加速壽命試驗(yàn)，應(yīng)力水平一般不少于4 個(gè)，每個(gè)應(yīng)力水平下的樣本量不少于5 個(gè)；Yang[5]提出了如何設(shè)計(jì)壽命滿足威布爾分布的4 個(gè)應(yīng)力水平加速壽命試驗(yàn)；一些學(xué)者基于蒙特卡羅方法提出了應(yīng)力水平、樣本量等設(shè)計(jì)要素的優(yōu)化設(shè)計(jì)方法[6-7]。針對(duì)加速壽命試驗(yàn)數(shù)據(jù)，主要采用統(tǒng)計(jì)學(xué)方法進(jìn)行分析，例如圖估計(jì)法[8]、簡單線性無偏估計(jì)法[9]、最好線性無偏估計(jì)法[10]、最小二乘法[11]、極大似然估計(jì)法[12]和貝葉斯法[13]。為了進(jìn)一步減小加速壽命試驗(yàn)的分析誤差，一些方法被用來對(duì)現(xiàn)有模型進(jìn)行修正。Mazzuchi 等[14]提出采用動(dòng)態(tài)線性模型與貝葉斯綜合方法來降低分析結(jié)果誤差；Watkins[15]提出一種修正的極大似然估計(jì)法，降低了計(jì)算的復(fù)雜度；殷毅超[16]提出了基于不精確概率理論的加速壽命試驗(yàn)數(shù)據(jù)分析方法。從工程應(yīng)用的角度分析，一些傳統(tǒng)的方法輔助工程經(jīng)驗(yàn)可以較好地解決單應(yīng)力多個(gè)應(yīng)力水平加速壽命試驗(yàn)設(shè)計(jì)和評(píng)估問題，而一些改進(jìn)和優(yōu)化方法仍停留在研究層面或者因求解方法過于復(fù)雜而難以推廣應(yīng)用。

在軟件工具方面，美國的ReliaSoft 軟件[17]、JMP Pro 軟件[18]、Minitab 軟件[19]和中國的ReliaQube軟件[20]均具有加速壽命試驗(yàn)?zāi)K，可以滿足單應(yīng)力多個(gè)應(yīng)力水平加速壽命試驗(yàn)設(shè)計(jì)和分析需求，但在適用性和易用性方面存在一些不足。以Minitab 軟件為例，它包含的加速模型僅有線性模型、自然對(duì)數(shù)模型等少數(shù)幾種，不包括阿倫尼斯模型等常用模型，并且在結(jié)果的可視化顯示方面存在不足。除上述軟件外，我國還有零星的用于研究目的的加速壽命試驗(yàn)軟件，但其研發(fā)和應(yīng)用尚處于起步階段，例如：重慶工業(yè)自動(dòng)化儀表研究所開發(fā)的加速壽命試驗(yàn)數(shù)據(jù)分析軟件[21]具備壽命分布擬合和加速模型參數(shù)估計(jì)等功能；電子科技大學(xué)開發(fā)了加速壽命試驗(yàn)方案優(yōu)化設(shè)計(jì)系統(tǒng)[22]，在試驗(yàn)條件設(shè)計(jì)方面很有特色，并融入了自主提出的改進(jìn)的隨機(jī)游走算法。然而，這些軟件在可視化顯示、工程經(jīng)驗(yàn)融合、分析結(jié)果驗(yàn)證等方面還存在一些不足，尚無實(shí)際的工程應(yīng)用案例。

本文將從軟件實(shí)現(xiàn)、計(jì)算方法和實(shí)例應(yīng)用的角度，對(duì)北京衛(wèi)星環(huán)境工程研究所自研的ALT511軟件[23]中單應(yīng)力多個(gè)應(yīng)力水平加速壽命試驗(yàn)設(shè)計(jì)與評(píng)估功能進(jìn)行系統(tǒng)和全面的闡述，以期對(duì)我國加速壽命試驗(yàn)軟件的設(shè)計(jì)和應(yīng)用提供有益參考和支撐。

1 單應(yīng)力多個(gè)應(yīng)力水平模塊介紹

1.1 “導(dǎo)入數(shù)據(jù)”選項(xiàng)卡

如圖1 所示，“導(dǎo)入數(shù)據(jù)”選項(xiàng)卡主要用于讀入各個(gè)應(yīng)力水平下的失效時(shí)間數(shù)據(jù)；點(diǎn)擊下方的“導(dǎo)入數(shù)據(jù)”按鈕會(huì)彈出選擇數(shù)據(jù)文件窗口，選定目標(biāo)數(shù)據(jù)文件后，數(shù)據(jù)將以表格形式進(jìn)行可視化呈現(xiàn)。本軟件可處理4 組應(yīng)力水平數(shù)據(jù)，并可兼容3 組和2 組應(yīng)力水平的情況。由于航天器產(chǎn)品通常價(jià)格昂貴且樣本量有限，所以多個(gè)應(yīng)力水平的加速壽命試驗(yàn)直接應(yīng)用于單機(jī)產(chǎn)品的情況并不多見，它更適用于涂層等材料以及軸承、電機(jī)等零件。應(yīng)力水平建議不少于4 個(gè)，每個(gè)應(yīng)力水平下的樣本量建議不少于5 個(gè)。

圖1 “導(dǎo)入數(shù)據(jù)” 選項(xiàng)卡界面Fig.1 Interface of the tab “data import”

1.2 “每組應(yīng)力水平壽命分布”選項(xiàng)卡

如圖2 所示，“每組應(yīng)力水平壽命分布”選項(xiàng)卡主要針對(duì)指數(shù)分布和威布爾分布進(jìn)行數(shù)據(jù)擬合并繪制概率密度函數(shù)；可針對(duì)指數(shù)分布給出各個(gè)應(yīng)力水平下的失效率，針對(duì)威布爾分布給出各個(gè)應(yīng)力水平下的形狀參數(shù)和特征壽命。

圖2 “每組應(yīng)力水平壽命分布”選項(xiàng)卡界面Fig.2 Interface of the tab “l(fā)ife distribution under each stress level”

分布類型通過界面左上方的下拉菜單選擇；根據(jù)選擇的分布類型，界面自動(dòng)更新顯示相應(yīng)的曲線和壽命分布模型參數(shù)。建議偶然性失效的產(chǎn)品（例如星載電源產(chǎn)品）選用指數(shù)分布；損耗性失效的產(chǎn)品（例如控制力矩陀螺軸承、艙外電動(dòng)工具按鈕、超聲電機(jī)、太陽帆板驅(qū)動(dòng)機(jī)構(gòu)導(dǎo)電滑環(huán)等）選用威布爾分布。

1.3 “加速模型擬合”選項(xiàng)卡

“加速模型擬合”選項(xiàng)卡可選擇逆冪律模型和阿倫尼斯模型2 種加速模型；針對(duì)指數(shù)分布可對(duì)不同應(yīng)力水平的平均故障前工作時(shí)間（MTTF）進(jìn)行擬合建模，針對(duì)威布爾分布可對(duì)不同應(yīng)力水平的特征壽命進(jìn)行擬合建模，最終將擬合結(jié)果繪制成曲線；繪圖界面中同時(shí)顯示擬合點(diǎn)和擬合線，同時(shí)根據(jù)最小二乘法計(jì)算并給出擬合相關(guān)系數(shù)，并顯示加速模型公式以及公式中模型常數(shù)的計(jì)算結(jié)果。需要指出，本軟件繪制加速模型曲線共4 種情況，即2 種壽命分布和2 種加速模型的排列組合，分別是指數(shù)分布逆冪律加速模型、指數(shù)分布阿倫尼斯加速模型、威布爾分布逆冪律加速模型和威布爾分布阿倫尼斯加速模型（如圖3 所示）。

圖3 “加速模型擬合”選項(xiàng)卡界面Fig.3 Interface of the tab “acceleration model fitting”

1.4 “指定應(yīng)力水平下壽命值”選項(xiàng)卡

“指定應(yīng)力水平下壽命值”選項(xiàng)卡主要用來計(jì)算可靠壽命、失效概率密度函數(shù)和可靠度函數(shù)。如圖4 所示，在界面中輸入指定應(yīng)力水平、可靠度、置信度，然后在下拉菜單中選擇是否考慮置信度，即可給出可靠壽命，并自動(dòng)繪制或更新對(duì)應(yīng)的概率密度函數(shù)曲線和可靠度函數(shù)曲線。

圖4 “指定應(yīng)力水平下壽命值”選項(xiàng)卡界面Fig.4 Interface of the tab “reliable life value under specified stress level”

2 計(jì)算方法

2.1 每組應(yīng)力水平壽命分布

1）指數(shù)分布

對(duì)于指數(shù)分布，主要進(jìn)行失效率參數(shù)估計(jì)，進(jìn)而得到失效概率密度函數(shù)[24]。每組應(yīng)力水平下失效率可根據(jù)極大似然估計(jì)法計(jì)算，

式中：下標(biāo)i指第i個(gè)應(yīng)力水平；Ti為累計(jì)試驗(yàn)時(shí)間；ri為失效數(shù)； λ?i為失效率估計(jì)值。平均失效前工作時(shí)間為

第i個(gè)應(yīng)力水平下失效概率密度函數(shù)為

2）威布爾分布

對(duì)于威布爾分布，主要進(jìn)行形狀參數(shù)和特征壽命兩個(gè)參數(shù)的估計(jì)，進(jìn)而得到失效概率密度函數(shù)[24]。對(duì)于失效數(shù)≥3 的情況，按照極大似然估計(jì)法通過求解方程(4)計(jì)算每組應(yīng)力水平下的形狀參數(shù)。

式中：mi為形狀參數(shù)；ri為失效數(shù)；tij為第j個(gè)試驗(yàn)件的試驗(yàn)時(shí)間；ni為樣本量。

對(duì)于失效數(shù)＜3 的情況，根據(jù)工程經(jīng)驗(yàn)確定形狀參數(shù)mi。第i個(gè)應(yīng)力水平下特征壽命估計(jì)值 η?i按照式(5)和式(6)計(jì)算：

第i個(gè)應(yīng)力水平下失效概率密度函數(shù)為

2.2 加速模型擬合

以下分別針對(duì)阿倫尼斯模型和逆冪律模型介紹如何擬合不同應(yīng)力水平下的壽命。

1）阿倫尼斯模型

用阿倫尼斯模型表達(dá)壽命與應(yīng)力之間的關(guān)系：

式中：a和b為模型常數(shù)；S為應(yīng)力；L為應(yīng)力水平是S時(shí)的壽命。這里的壽命L對(duì)于指數(shù)分布為平均失效前工作時(shí)間θ，對(duì)于威布爾分布為特征壽命η。對(duì)式(8)等號(hào)兩邊取對(duì)數(shù)進(jìn)行線性化處理得

式中：Y=lnL；A=lna；B=b；X=1/S。根據(jù)最小二乘法，可得常數(shù)B和A的估計(jì)值分別為：

其中：Xi和Yi分別為第i個(gè)應(yīng)力水平下的X和Y取值；Xˉ 和Yˉ分別為X和Y的平均值；k為應(yīng)力水平的數(shù)量。根據(jù)計(jì)算出來的A? 和B?可以反求加速模型的常數(shù)a和b。擬合相關(guān)系數(shù)ρ為

其中，ρ越接近1，代表選擇該加速模型與試驗(yàn)數(shù)據(jù)的匹配程度越高。

2）逆冪律模型

用逆冪律模型表達(dá)的壽命與應(yīng)力之間的關(guān)系為

式中：c和α為模型常數(shù)。同樣，這里的壽命對(duì)于指數(shù)分布為平均失效前工作時(shí)間 θ?i，對(duì)于威布爾分布為特征壽命 η?i。對(duì)式(13)等號(hào)兩邊取對(duì)數(shù)同樣可得到形如式(9)的表達(dá)式，其中，Y=lnL，A=lnc，B=-α，X=lnS。類似阿倫尼斯模型擬合，可按照式(10)和式(11)分別計(jì)算常數(shù)B和A，并反求加速模型的常數(shù)c和α。進(jìn)而可根據(jù)式(12)計(jì)算擬合相關(guān)系數(shù)，評(píng)估逆冪律模型與試驗(yàn)數(shù)據(jù)的匹配程度。

2.3 指定應(yīng)力水平下壽命值

根據(jù)2.2 節(jié)得到的加速模型，可計(jì)算指定應(yīng)力水平下的壽命。這里的壽命對(duì)于指數(shù)分布為平均失效前工作時(shí)間θ，對(duì)于威布爾分布為特征壽命η。以下介紹在不考慮置信度和考慮置信度下，如何得到失效概率密度函數(shù)、可靠度函數(shù)和可靠壽命。

1）不考慮置信度

在不考慮置信度的情況下，可以得到指數(shù)分布的失效概率密度函數(shù)f(t)和可靠度函數(shù)R(t)表達(dá)式，分別見式(14)和式(15)，并可繪制曲線。

威布爾分布的失效概率密度函數(shù)f(t)和可靠度函數(shù)R(t)表達(dá)式分別見式(16)和式(17)，并可繪制曲線。

指數(shù)分布和威布爾分布的可靠壽命tR表達(dá)式可分別由式(15)和式(17)推導(dǎo)得出，分別為

2）考慮置信度

在考慮置信度的情況下，對(duì)于指數(shù)分布，采用式(20)[24]得到平均失效前工作時(shí)間的置信下限θL：

圖5 置信壽命系數(shù)與失效數(shù)、置信度、形狀參數(shù)的關(guān)系Fig.5 Impact of failure numbers, confidence level and shape parameter on confidence life factor

考慮置信度的失效概率上限密度函數(shù)fu(t)和可靠度下限函數(shù)RL(t)表達(dá)式分別為式(21)和式(22)，可靠壽命下限tRL表達(dá)式為式(23)。

對(duì)于威布爾分布，采用式(24)得到特征壽命的置信下限ηL：

對(duì)比式(24)和式(20)，可以發(fā)現(xiàn)指數(shù)分布是威布爾分布形狀參數(shù)mˉ取1 的特殊情況。將式(24)中ηL/η記為威布爾分布的置信壽命系數(shù)，表征考慮置信度下特征壽命的降低程度，從圖5(b)～(d)可以看出：置信壽命系數(shù)隨著失效數(shù)的增加而增大、隨著置信度的增加而減小，這一規(guī)律與指數(shù)分布的一致；并且隨著mˉ增加，置信壽命系數(shù)顯著降低。

類似指數(shù)分布，可以得到考慮置信度的失效概率上限密度函數(shù)fu(t)和可靠度下限函數(shù)RL(t)表達(dá)式，分別見式(25)和式(26)；可靠壽命下限tRL表達(dá)式見式(27)。

3 軟件應(yīng)用案例

以某航天器超聲電機(jī)為對(duì)象，給出單應(yīng)力多個(gè)應(yīng)力水平加速壽命試驗(yàn)設(shè)計(jì)和評(píng)估的軟件應(yīng)用案例。

根據(jù)失效模式與影響分析（FMEA），以負(fù)載為加速應(yīng)力，設(shè)計(jì)比正常應(yīng)力水平0.02 N·m 高且比破壞極限0.22 N·m 低的4 個(gè)應(yīng)力水平開展加速壽命試驗(yàn)并記錄試驗(yàn)件的失效時(shí)間。如圖6 所示，將0.03 、0.07 、0.16 和0.21 N·m 共4 個(gè)應(yīng)力水平下共20 個(gè)試驗(yàn)件的失效數(shù)據(jù)讀入軟件并顯示。

圖6 某航天器超聲電機(jī)在4 個(gè)應(yīng)力水平下的失效時(shí)間Fig.6 Failure time of a spacecraft ultrasonic motor under four stress levels

又根據(jù)FMEA，該航天器超聲電機(jī)在負(fù)載應(yīng)力作用下主要失效機(jī)理為磨損失效，具有典型的累積損傷失效特征，這樣各個(gè)應(yīng)力水平下的壽命應(yīng)滿足威布爾分布。選擇分布類型后，軟件可自動(dòng)計(jì)算4 個(gè)應(yīng)力水平下的形狀參數(shù)、特征壽命，并繪制失效概率密度函數(shù)曲線（如圖7 所示）。以應(yīng)力水平1（0.03 N·m）為例，壽命滿足形狀參數(shù)為5.1、特征壽命為5941 h 的威布爾分布。此外，觀察4 個(gè)應(yīng)力水平的壽命分布，隨著負(fù)載的增大，特征壽命逐漸降低，說明加大負(fù)載能很好起到縮短試驗(yàn)時(shí)間的加速效果；而4 個(gè)應(yīng)力水平的形狀參數(shù)相差不大，提示各個(gè)應(yīng)力水平的失效機(jī)理基本一致，各組應(yīng)力水平下的試驗(yàn)數(shù)據(jù)有效，可用于建立加速模型。

圖7 4 組應(yīng)力水平下的某航天器超聲電機(jī)壽命分布Fig.7 Life distribution of a spacecraft ultrasonic motor under four stress levels

基于4 個(gè)應(yīng)力水平下的特征壽命數(shù)據(jù)進(jìn)行加速模型擬合，結(jié)果如圖8 所示，可見：當(dāng)選擇阿倫尼斯模型時(shí)，算得模型常數(shù)a和b分別為32.22 和0.171，相關(guān)系數(shù)為0.894 4；而選擇逆冪律模型時(shí)，算得模型常數(shù)c和α分別為0.519 3 和2.776，相關(guān)系數(shù)為0.971 9。由于逆冪律模型的相關(guān)系數(shù)更接近1，說明數(shù)據(jù)擬合的效果更好，所以選擇逆冪律模型作為航天器超聲電機(jī)的加速模型。這一加速模型表達(dá)式與圖9 所示Excel 軟件的擬合結(jié)果一致，證明加速模型擬合算法是正確的。但Excel 軟件沒有針對(duì)阿倫尼斯模型的擬合功能，Minitab 軟件的加速壽命試驗(yàn)?zāi)K也未提供阿倫尼斯加速模型和逆冪律加速模型的擬合功能，這從一個(gè)側(cè)面反襯出本軟件的專業(yè)性。

圖8 某航天器超聲電機(jī)加速模型擬合Fig.8 Acceleration model fitting for a spacecraft ultrasonic motor

圖9 Excel 軟件擬合的逆冪律加速模型Fig.9 Inverse power law acceleration model fitted by Excel

獲得加速模型后，可得到指定應(yīng)力水平、可靠度和置信度要求下的可靠壽命。如圖10 所示，假設(shè)0.02 N·m 為正常應(yīng)力水平，這一負(fù)載下，不考慮置信度時(shí)可靠度0.95 對(duì)應(yīng)的可靠壽命為15 180 h，這一計(jì)算結(jié)果與圖11 所示的Minitab 軟件計(jì)算結(jié)果幾乎一致。

圖10 某航天器超聲電機(jī)可靠壽命計(jì)算Fig.10 Calculation of the reliable life of a spacecraft ultrasonic motor

圖11 Minitab 軟件計(jì)算的可靠度0.95 的可靠壽命Fig.11 Reliable life under 0.95 reliability calculated by Minitab

此外，本軟件提供了更豐富的概率密度函數(shù)曲線和可靠度函數(shù)曲線顯示功能，例如：設(shè)置置信度0.7，考慮置信度條件下，本軟件計(jì)算的可靠度0.95 對(duì)應(yīng)的可靠壽命為14 760 h。加速壽命試驗(yàn)中各試驗(yàn)件最長試驗(yàn)時(shí)長為7100 h、最短試驗(yàn)時(shí)長12.5 h，均遠(yuǎn)小于要驗(yàn)證的壽命指標(biāo)要求。而且完成多個(gè)應(yīng)力水平的加速壽命試驗(yàn)后，后續(xù)針對(duì)同類產(chǎn)品可以只做一個(gè)短耗時(shí)的高應(yīng)力水平加速壽命試驗(yàn)，即可驗(yàn)證壽命指標(biāo)。

4 結(jié)束語

本文介紹了加速壽命試驗(yàn)設(shè)計(jì)與評(píng)估軟件ALT511 的單應(yīng)力多個(gè)應(yīng)力水平模塊，包括：詳細(xì)說明了“導(dǎo)入數(shù)據(jù)”、“每組應(yīng)力水平壽命分布”、“加速模型擬合”和“指定應(yīng)力水平下壽命值”4 個(gè)選項(xiàng)卡；給出了基于多應(yīng)力水平失效數(shù)據(jù)的指數(shù)分布和威布爾分布擬合方法、阿倫尼斯模型和逆冪律模型的加速模型擬合方法，以及指定應(yīng)力水平下考慮和不考慮置信度兩種情況的可靠壽命、概率密度函數(shù)和可靠度函數(shù)計(jì)算方法。以某航天器超聲電機(jī)加速壽命試驗(yàn)為例，軟件應(yīng)用結(jié)果表明：ALT511 軟件的加速模型擬合結(jié)果與Excel 軟件擬合結(jié)果一致，逆冪律模型常數(shù)c和α分別為0.519 3 和2.776，相關(guān)系數(shù)為0.971 9；可靠壽命計(jì)算結(jié)果與Minitab 軟件計(jì)算結(jié)果幾乎一致。

后續(xù)還將進(jìn)一步研究多應(yīng)力加速壽命試驗(yàn)、加速常數(shù)靈敏度分析、加速模型可信度分析，以及這些技術(shù)的軟件實(shí)現(xiàn)。