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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

1.3 “加速模型擬合”選項卡

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

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

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

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

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

2 計算方法

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

1）指數(shù)分布

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

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

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

2）威布爾分布

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

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

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

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

2.2 加速模型擬合

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

1）阿倫尼斯模型

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

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

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

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

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

2）逆冪律模型

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

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

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

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

在考慮置信度的情況下，對于指數(shù)分布，采用式(20)[24]得到平均失效前工作時間的置信下限θ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)。

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

對比式(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)用案例

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

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

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

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

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

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

圖8 某航天器超聲電機加速模型擬合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ù)載下，不考慮置信度時可靠度0.95 對應(yīng)的可靠壽命為15 180 h，這一計算結(jié)果與圖11 所示的Minitab 軟件計算結(jié)果幾乎一致。

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

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

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

4 結(jié)束語

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

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