宋征宇,方志耕,赫武樂,孫云柯,王 召,李彩霞,劉思峰
(1. 中國運載火箭技術研究院,北京 100076;2. 南京航空航天大學 經(jīng)濟與管理學院,南京 211106)
運載火箭研制具有高成本、高可靠和高風險等特點。為保證執(zhí)行任務足夠高的可靠性,開展可靠性分析、評估等量化工作尤為關鍵。選取恰當?shù)目煽啃栽u估方法是運載火箭研制工程的迫切需求。
經(jīng)典可靠性評估理論和方法是基于概率論和數(shù)理統(tǒng)計學建立的,由產(chǎn)品可靠性試驗數(shù)據(jù)統(tǒng)計推斷產(chǎn)品的可靠性。基于概率論的可靠性評估,需要首先擬合產(chǎn)品壽命分布。從實際來看,電子產(chǎn)品、機械產(chǎn)品和機電產(chǎn)品等不同類型產(chǎn)品服從不同的壽命分布模型。針對常見單元分布類型的經(jīng)典可靠性評估方法,Epstein B.等,Weibull W.,Daniels H.E.,分別對指數(shù)型、Weibull型、正態(tài)型、對數(shù)正態(tài)型單元的評估方法進行了研究[1-3]?,F(xiàn)有運載火箭單機產(chǎn)品可靠性評估模型主要包括指數(shù)壽命、性能正態(tài)等評估方法、應力–強度法等。呂箴等在對運載火箭典型產(chǎn)品常用可靠性評估特征量提取的基礎上,針對新一代運載火箭機構類產(chǎn)品構成復雜、失效模式種類繁多,部分產(chǎn)品評估信息少的問題,提出“三要素”可靠性評估方法[4]。申桂香等[5]采用極大似然估計與參數(shù)偏差修正方法,對威布爾分布中的參數(shù)進行了估計與修正,對于數(shù)控機床可靠性評估中的小樣本類型故障數(shù)據(jù)處理方法進行了研究,結果對可靠性評估模型的偏差修正具有顯著效果。平仕良等[6]結合結構動態(tài)應力、振動等量化參數(shù),在考慮運載火箭發(fā)射臺的熱燒蝕、振動環(huán)境等影響下,采用專家模糊評分和層次分析法,對運載火箭發(fā)射臺主體結構可靠性進行了評估。馮鐵山等[7]以試驗階段和相似系統(tǒng)飛行試驗的有關歷史數(shù)據(jù)為先驗信息,綜合現(xiàn)場試驗信息,運用貝葉斯公式對運載火箭系統(tǒng)可靠性進行了評估。根據(jù)已有對運載火箭可靠性評估的研究可知,在對運載火箭首飛可靠性評估的過程中,受限于試驗周期短,試驗樣本成本高等因素,可靠性試驗子樣相對較少,數(shù)據(jù)較為匱乏,因此運載火箭的可靠性評估是典型的小子樣可靠性評估問題。
在國內(nèi)外研究中,小子樣可靠性評定方法在不斷涌現(xiàn),有置信區(qū)間法、貝葉斯方法、信息熵法、近似正態(tài)法、矩估計法等[8-11]。Hathout[12]將概率論和模糊理論相結合,應用于可靠性和安全性評估問題;Bourinet等將信息熵法應用于離散狀態(tài)空間馬爾可夫鏈模型,用交叉熵逐步更新傳遞概率矩陣,改進了蒙特卡羅方法在可靠性評估中的應用[13];文獻[14]提出了一種利用少量有價值產(chǎn)品的試驗數(shù)據(jù)進行可靠性評估的方法,將仿真或數(shù)字設計結果與實驗數(shù)據(jù)相結合,推斷出多狀態(tài)(多種群)下的實驗數(shù)據(jù),該方法不僅可以估計每個狀態(tài)下的均值和方差,還可以估計總體百分位數(shù)和百分比的置信限和區(qū)間,與傳統(tǒng)的可靠性評估方法相比,該方法具有較高的精度,解決了單一狀態(tài)下的可靠性評估問題;文獻[15]針對核工程和航天工程中的小樣本高可靠性安全失效模式器件的可靠性評估難題,分析了貝葉斯統(tǒng)計、改進的貝葉斯方法和貝葉斯網(wǎng)絡等可靠性評估方法,出于計算量、建模及信息處理方面的考慮,在實際可靠性工程中常用貝葉斯方法來解決小子樣復雜裝備可靠性評估,它是是解決可靠性信息不足的有效方法,能夠較好地綜合利用各種主、客觀的先驗信息以及多層的試驗數(shù)據(jù),相比傳統(tǒng)統(tǒng)計方法減少了所需試驗次數(shù)、縮短了試驗周期,節(jié)約了試驗成本,提高了試驗效率。美軍在80年代就采用貝葉斯小子樣理論進行“潘興Ⅱ”導彈的可靠性評定。1984年,美軍提出:對于昂貴武器系統(tǒng)進行破壞性試驗必須采用序貫分析方法或貝葉斯小子樣理論。俄羅斯和法國在進行系統(tǒng)可靠性評定時,也十分強調(diào)充分利用補充信息的貝葉斯方法或貝葉斯經(jīng)驗方法[16]。綜合上述研究可以發(fā)現(xiàn),現(xiàn)有方法大多以工程經(jīng)驗信息和歷史試驗信息作為先驗信息,并結合現(xiàn)場試驗信息對產(chǎn)品可靠性進行評估。
貝葉斯方法能夠較好地綜合利用各種主、客觀的先驗信息以及多層的試驗數(shù)據(jù),解決可靠性信息不足的問題,在實際可靠性工程中常被用于小子樣復雜裝備的可靠性評估。因此,近年來小子樣系統(tǒng)的Bayes可靠性綜合方法在國內(nèi)外得到了更多的研究和應用[17-20]。郭凱紅[21]介紹了基于貝葉斯理論的小子樣可靠性評估的流程、數(shù)學模型、數(shù)值計算方法,建立了可靠性評估應用方案,有效指導了該系統(tǒng)的可靠性評估工作,為其它型號開展小子樣可靠性評估應用提供借鑒。羅潤[22]結合工程經(jīng)驗,介紹了貝葉斯方法在飛機可靠性評估工作中的應用。李婧等[23]將貝葉斯理論運用于武器裝備的貯存可靠性評估問題,結合貯存期間檢測數(shù)據(jù),構建了基于貝葉斯推斷的指數(shù)型單元貯存效果評估模型,對指數(shù)型單元貯存效果實現(xiàn)進行了合理的評估。周源泉等[24]進行了“長征”系列運載火箭可靠性增長分析,并分別用經(jīng)典和貝葉斯方法對當前LM火箭的可靠性進行了評估,驗證了貝葉斯方法在運載火箭可靠性評估方面的有效性。伯仲干等[25]針對不同分布類型的單元,對成敗型、指數(shù)型和weibull型分布單元貝葉斯可靠性評估方法進行了研究,并提出了基于信度加權的單元可靠性數(shù)據(jù)融合評估新方法,對傳統(tǒng)的貝葉斯可靠性評估方法作出了改進。李大偉等[26]針對現(xiàn)場使用可靠性信息收集困難的問題,利用貝葉斯方法處理產(chǎn)品備件需求信息,對長期貯存產(chǎn)品的使用可靠性進行評估。邵松世等[27]通過定義不同先驗信息的似然函數(shù)系數(shù),進行了貝葉斯可靠性評估中多源先驗數(shù)據(jù)的融合方法研究,建立了指數(shù)型備件可靠性的貝葉斯評估方法,提升了評估結果的準確度與穩(wěn)定性。綜合上述研究,近年來小子樣系統(tǒng)的貝葉斯可靠性綜合方法得到了很多研究和應用,特別是在有較多的歷史數(shù)據(jù)或較強的主觀信息的場合,不僅可以節(jié)省大量的試驗經(jīng)費,還可縮短整個系統(tǒng)的研制周期,已成為可靠性工程師和統(tǒng)計工作者普遍研究的方法。但現(xiàn)有的小子樣貝葉斯可靠性評估方法的共軛先驗分布獲取基本還建立在經(jīng)典的共軛先驗分布表基礎之上,針對威布爾分布這類總體分布,由于沒有共軛先驗分布,先驗分布的獲取極為困難,貝葉斯推導繁雜且大概率無解。
本文針對威布爾、二項、指數(shù)3個可靠性評估工作中的典型總體分布,采取新的建模策略,分別構建相應的可靠性評估貝葉斯模型,解決威布爾分布沒有共軛先驗分布的總體分布貝葉斯建模與求解問題。本文充分利用歷史統(tǒng)計數(shù)據(jù)、經(jīng)驗信息及小樣本數(shù)據(jù),做出設備可靠性參數(shù)的貝葉斯估計。通過數(shù)據(jù)處理、模型構建、求解等過程,總結了基于貝葉斯理論的小子樣可靠性評估流程。最后,對運載火箭相關設備可靠度進行評估,通過案例分析、方法對比,驗證本文方法的有效性與優(yōu)越性。
貝葉斯公式的一般形式為
其中:p(x|θ)為隨機變量x的似然函數(shù);π(θ)為參數(shù)θ 的先驗分布密度;π(θ|x)為參數(shù)θ 的后驗分布密度。
當θ 是離散隨機變量時,貝葉斯公式為
由貝葉斯公式可知:在試驗X=x前,對參數(shù)的認識總結于π(θ)中,而試驗X=x所取得的關于的新信息則包含在似然函數(shù)p(x|θ)中,經(jīng)修正后,先驗分布π(θ)變?yōu)楹篁灧植鸡?(θ|x),即經(jīng)實踐后,修正了原有的認識,達到了更高一級的認識。貝葉斯公式中后驗分布是綜合了先驗分布信息和樣本信息。
事實證明,定理的分母項不依賴于參數(shù)θ ,因此貝葉斯公式也可以表示為
威布爾分布對各種類型試驗數(shù)據(jù)擬合能力極強,能擬合運載火箭中大多數(shù)設備、系統(tǒng)的壽命試驗數(shù)據(jù),構建二參數(shù)威布爾壽命分布設備的貝葉斯可靠性評估模型為
1)威布爾總體分布
總體分布為威布爾分布,其分布函數(shù)為
其中:m為形狀參數(shù);η 為設備的特征壽命。
可靠度函數(shù)為
令λ=η?m,則式(5)等價于
由于在可靠性評定過程中,威布爾參數(shù)的共軛先驗分布不存在,因此考慮其可靠性函數(shù)為指數(shù)函數(shù),采用指數(shù)分布的共軛先驗分布,即伽馬分布擬合先驗數(shù)據(jù)信息,且伽馬分布的分布族非常豐富,基本能擬合絕大部分分布情況,其表達式易于進行貝葉斯推斷。因此將威布爾分布轉化為關于統(tǒng)計量 λ的指數(shù)函數(shù),由表1可知,指數(shù)函數(shù)的共軛先驗分布為伽馬分布,則取λ 的共軛先驗分布進行下一步的計算。
表1 常見的共軛先驗分布Table 1 The usual conjugate prior distribution
假設
根據(jù)設備可靠性先驗信息,以矩估計的方法求解(a,b),解方程組(8)
即有統(tǒng)計量λ 的先驗分布
3)總體、樣本信息的似然函數(shù)
n個產(chǎn)品進行試驗,截尾時間為t,則壽命服從威布爾在無失效的情況下,λ 的似然函數(shù)為
4)統(tǒng)計量λ 的后驗分布
根據(jù)貝葉斯公式結合公式(8)、(9),可得λ 的后驗密度為
5)設備特征壽命的置信下限計算
設特征壽命η 的置信下限為ηL,即滿足
則一定任務時間t的可靠度為
其中:m為形狀參數(shù),取2或2.4; η為特征壽命;t為規(guī)定任務時間。
成敗型單元可靠度總體分布服從二項分布,參見表1,工程上通常以 β分布作為共軛先驗分布,求解后驗分布,評估產(chǎn)品可靠度。
1)先驗分布:f0(R)=β(R|s0,f0),由先驗數(shù)據(jù)可得成功次數(shù)s0和失敗次數(shù)f0。
3)后驗密度則由貝葉斯定理確定
4)在置信度 γ的可靠度置信下限RL,B為
由此可求得可靠度置信下限RL,B。
指數(shù)分布描述了產(chǎn)品處在隨機失效階段時的壽命分布,由于這個階段占據(jù)了產(chǎn)品的大部分工作時間,因此指數(shù)壽命型分布是可靠性工程領域最常見的一種分布類型。
在指數(shù)型設備的可靠性評估中,對于設備的檢驗方式有定時截尾(有替換、無替換),定數(shù)截尾(有替換、無替換)以及隨機截尾的檢驗方式。主要評估的指標有故障率 λ的上限λu或平均壽命θ (θ = λ?1)的下限 θL,可靠度R(R=e?λt0)下限RL以及給定可靠度的壽命tR(tR= λ?1ln(1/R))下限tR,L,以下利用有替換定時截尾試驗數(shù)據(jù),構建指數(shù)型部件可靠度貝葉斯模型[28]。
1)指數(shù)型總體分布
從壽命服從指數(shù)分布F(t)=1?e?λt的母體中抽樣,作n個產(chǎn)品的有替換定時截尾試驗,總試驗時間τ,試驗時間內(nèi)失效數(shù)為z,規(guī)定任務時間為t,其失效率服從分布。
2)失效率λ、可靠度R的先驗分布
由表1可知,取失效率λ共軛先驗分布為伽馬分布,其先驗密度函數(shù)為
其中:z0為形狀參數(shù);τ0為尺度參數(shù)。
則根據(jù)變換下的不變性原則可得該單元的任務時間的可靠度R=R(t)=e?λt的先驗分布為
即R的共軛驗前分布為負對數(shù)伽瑪分布,記為LΓ(z0,τ0/t),并記 β =τ0/t。
3)先驗超參數(shù)的確定
z0和τ0可由先驗信息直接確定,即在總時間為τ0的試驗中失效數(shù)為z0。也可在已知參數(shù)λ的驗前期望與方差的情況下,通過求解方程組(19)來確定
4)可靠度R的后驗分布
根據(jù)貝葉斯定理,代入現(xiàn)場試驗樣本,整理可得到后驗分布如下
5)設備可靠度的置信下限計算
將 π (R|D)代入下面式子,即可求出產(chǎn)品在給定置信度γ下的置信下限RL
以貝葉斯方法為主,經(jīng)典方法為輔,通過前期對試驗數(shù)據(jù)的收集與預處理,結合不同類型設備的工作原理,利用貝葉斯定理可最終得出產(chǎn)品(包括系統(tǒng)、分系統(tǒng)、單元)的可靠性、貝葉斯可靠性下限等參數(shù)的評估結果,驗證產(chǎn)品是否達到了分配和預計的可靠性指標,形成一套完整的產(chǎn)品可靠性評估方案。
對產(chǎn)品(單元、分系統(tǒng)、系統(tǒng))的可靠性數(shù)據(jù)進行完善的記錄并進行初步的統(tǒng)計分析??煽啃詳?shù)據(jù)包括試驗數(shù)據(jù)和現(xiàn)場數(shù)據(jù)。試驗數(shù)據(jù)來自可靠性試驗、壽命試驗、加速壽命試驗、功能試驗、環(huán)境試驗和綜合試驗等??煽啃栽囼炛饕越匚苍囼灋橹鳎话銥槎〝?shù)截尾試驗、定時截尾試驗和隨機截尾試驗等。運載火箭設備單元主要有電子產(chǎn)品、機械產(chǎn)品和機電產(chǎn)品,單元設備的可靠性分布模型一般有兩項、指數(shù)、Weibull、正態(tài)和對數(shù)正態(tài)模型。
先驗信息主要來源于單元測試的歷史信息、類似單機產(chǎn)品的測試信息、單機仿真數(shù)據(jù)以及專家經(jīng)驗信息等,以上信息具有不確定性和一定程度上的主觀性,為使提出的方法具有嚴密的理論基礎,應首先建立先驗信息的統(tǒng)計信息。
可靠性信息可以分布函數(shù)的形式體現(xiàn),這是較為理想的表達方式。對于貝葉斯理論在可靠性評估中的應用,合理選取先驗分布才能保證評估結果有足夠穩(wěn)健性和可信度。構建先驗分布的方法有多種,這里主要應用最方便、廣泛的共軛先驗分布。
共軛先驗分布概念清晰,相關參數(shù)可由相對最小二乘或歷史試驗數(shù)據(jù)確定,計算簡便。因此共軛先驗分布也是工程上應用最廣泛的。表1給出了常見共軛先驗分布。共軛先驗分布的提出解決了貝葉斯推斷中后驗分布難以推出的難點,可以通過簡單的數(shù)字計算得到后驗分布。同時,常見的共軛先驗分布,分布族豐富,能擬合絕大多數(shù)先驗信息。
其它常見的方法還有最大熵先驗分布、蒙特卡羅法、Gibbs抽樣法、Bootstrap及隨機加權法。
以上4種先驗分布確定方法中,共軛先驗分布和最大熵方法適用于結構較為簡單的單元或系統(tǒng),在許多研究中將這兩種方法結合使用,利用最大熵方法求解共軛先驗分布中的超參數(shù),既可保持共軛先驗分布的“繼承性”優(yōu)勢,又可減少主觀因素的滲入,使計算結果更加精確;蒙特卡羅、Gibbs抽樣法及Bootstrap隨機加權法適用于結構較為復雜的單元或系統(tǒng),Bootstrap與隨機加權法也適用于小樣本情況下的先驗分布構建。
步驟3:后驗分布的確定
結合貝葉斯定理,參數(shù)的后驗分布正比于其先驗分布與試驗結果的似然函數(shù)之積,似然函數(shù)根據(jù)先驗信息較容易求得。根據(jù)后驗分布,可求得設備的可靠度、貝葉斯可靠性下限等參數(shù),得到指標驗證結論。
利用貝葉斯方法的設備可靠性評定流程如圖1所示:
圖1 基于貝葉斯理論的設備可靠性評定流程Fig. 1 The process of equipment reliability evaluation based on the Bayesian theory
以某型號運載火箭為例,該系統(tǒng)包括電子產(chǎn)品、機械產(chǎn)品和機電產(chǎn)品等服從不同壽命分布模型的單元。本節(jié)選取3種典型的運載火箭設備,結合上文提出的基于貝葉斯理論的可靠性評估方法,根據(jù)歷史經(jīng)驗以及試驗數(shù)據(jù)等先驗信息,構建不同工作環(huán)境下相應分布類型的可靠性評估模型,以對其進行可靠性評定。
3.1.1 飛行可靠度
1)試驗數(shù)據(jù)及已知信息
氧地面增壓管屬于機械產(chǎn)品,其壽命服從威布爾分布,現(xiàn)對氧地面增壓管進行可靠性評定。經(jīng)過關聯(lián)聚類分析,選取類似產(chǎn)品氫地面增壓管(試驗時間18 540 s),將類氫地面增壓管試驗信息(試驗時間18 540 s)作為先驗信息,代入威布爾分布可靠度計算公式,計算得先驗可靠度,由可靠性基本理論,可計算其特征壽命ηL(可靠度為0.368時對應的壽命)以及壽命方差σ(η)2,步驟如下
2)給出先驗分布
從先驗信息可知氧地面增壓管的特征壽命為100,壽命方差約為10,代入公式(8)
可得a= 0.000 000 1,b= 0.000 011
3)計算可靠度置信下限
在置信度為0.7的情況下,將已知數(shù)據(jù)a、b代入式(10),計算設備特征壽命的置信度下限
4)規(guī)定任務時間內(nèi)可靠度計算
在規(guī)定任務時間為1 min的條件下,將特征壽命置信下限代入式(14)求解設備飛行可靠度
解得Rt=0.995994790960113。
3.2.1 發(fā)射可靠度
A100發(fā)動機屬于機電產(chǎn)品,其點火試驗數(shù)據(jù)服從二項分布,以貝塔分布作為先驗分布,求解后驗分布,評估產(chǎn)品發(fā)射可靠度。
1)試驗數(shù)據(jù)及已知信息
現(xiàn)對A100發(fā)動機進行可靠性評定,已檢驗其點火成功與否服從二項分布。以A115發(fā)動機作為先驗試驗數(shù)據(jù),其總共實驗40次,失效次數(shù)f0= 1,成功次數(shù)s0= 39。
2)給出先驗分布
以貝塔分布作為先驗分布,將已知數(shù)據(jù)代入公式求得共軛型先驗密度函數(shù)
其發(fā)生的概率為
3)確定后驗密度
A100發(fā)動機點火成功次數(shù)284次,失效0次,其二項實驗結果(s,f)=(284,0),則A100發(fā)動機的后驗密度由貝葉斯定理確定
4)可靠度置信下限計算
置信度為0.7的可靠度置信下限為
由Excel中Beta.inv函數(shù)計算結果為
碳化硅功率器件屬于航天電子產(chǎn)品,已檢驗其壽命服從指數(shù)分布,對其進行定時截尾試驗,以負對數(shù)伽瑪分布作為先驗分布,求解后驗分布,評估該部件可靠度。
1)試驗數(shù)據(jù)及已知信息
指數(shù)單元作有替換定時截尾試驗,失敗數(shù)z=1,總試驗時間τ=100 h,任務工作時間t=10 h,根據(jù)歷史試驗信息可知z0=2,τ0=200h。
2)計算先驗分布
將已知數(shù)據(jù)代入公式(18)求得共軛型先驗密度函數(shù)
3)計算后驗密度
根據(jù)貝葉斯定理,將現(xiàn)場試驗樣本數(shù)據(jù)代入公式(20),整理可得到后驗分布
4)計算給定置信度下的可靠度下限
根據(jù)后驗密度及現(xiàn)場試驗信息,由公式(21)可計算該部件置信度0.9下的可靠度置信下限為
計算可得RL= 0.837 4。
表2給出了在相同置信度下,不同方法對于以上3種產(chǎn)品的可靠度評估結果。由表2可知,最小二乘估計結果低于仿真方法與Bayes方法,不符合運載火箭相關產(chǎn)品高可靠度要求的實際情況,且小子樣數(shù)據(jù)的情況不適用于線性回歸的方法。仿真方法與本文計算結果較為接近,但試驗時間較長,復雜度高,綜合以上結果,Bayes評估方法可以充分利用先驗信息,對樣本量要求較低,適用于小子樣運載火箭可靠性評估,且能夠節(jié)約試驗時間與計算成本。
表2 不同方法可靠度計算結果Table 2 Results of different methods
通過對上述3類設備的可靠度計算結果,可驗證貝葉斯方法能夠充分利用先驗信息,節(jié)省時間和經(jīng)費,而且分析方法程式化,易于工程人員掌握,計算結果顯示:基于貝葉斯方法的可靠性評定在引入先驗信息后,能夠有效提升小子樣產(chǎn)品的可靠性,解決了小子樣產(chǎn)品以經(jīng)典方法評估可靠性時,評估結果過低、與實際認知明顯不符、評估不準確等問題。因此貝葉斯方法明顯優(yōu)于經(jīng)典可靠性評估方法,在數(shù)據(jù)較少的小子樣可靠性評估問題中應用效果較好。該方法有利于提高單元設備可靠性評估的準確度與穩(wěn)定度,為運載火箭系統(tǒng)可靠性評價提供更為準確的數(shù)據(jù)基礎。
在工程實際中,有時會將試驗條件變化微小、試驗樣品不變的先驗數(shù)據(jù)歸納到樣本數(shù)據(jù)中以充實樣本數(shù)據(jù)量,這樣得到的可靠性評估結果將會有所增大(特定分布下可能不變,如二項分布)。因此,在選擇樣本數(shù)據(jù)和先驗數(shù)據(jù)時,應充分考慮分布的類型,對于如威布爾分布這類先驗數(shù)據(jù)與樣本數(shù)據(jù)劃分不同而對評估結果有影響的總體分布,應該充分將與樣本信息的試驗條件、試驗樣品變化微小的先驗數(shù)據(jù)歸為樣本數(shù)據(jù),以提升評估精度。下一步將考慮多源先驗信息的數(shù)據(jù)融合、分布函數(shù)中的參數(shù)確定與偏差修正、不同試驗條件下的置信度等問題,對該評價方法開展進一步的研究工作。