王博, 蔣平, 郭波
(國防科技大學 系統工程學院, 湖南 長沙 410073)
當前,國家對國防軍工裝備提出更高層次的要求,對重點關注的航天產品也加大了投入與重視。航天閥門作為液體火箭發動機的關鍵部件,要具備很高的可靠性水平,在研制階段就需準確評估其可靠性。航天閥門的工作原理與結構復雜,評估難度大,同時,由于小樣本、高可靠性的特點,在研制階段試驗中僅能得到極少失效甚至無失效的試驗數據。
航天閥門的失效模式和影響分析已經趨向成熟,針對閥門啟動(關閉)失靈、卡滯、閥門污染(泄漏)等故障模式的狀態識別技術、失效機理分析以及設計改進,已有大量研究成果。但是,在發動機的整機試車時,仍然會出現卡滯、泄漏等問題,這說明現有研究仍未充分考察實際工作環境,也難以準確評估閥門的可靠性。液體火箭發動機工作環境惡劣,閥門可靠性要求閥門在工況溫度和振動環境下,額定的時間內完成規定的開閉次數。因此航天閥門在定型之前要進行溫度、振動和動作等研制試驗,保證閥門在真實工作環境下的可靠性。在現有閥門可靠性評估成果中,這些研制試驗數據未得到充分利用。
本文提出融合多類型研制試驗數據的方法來評估航天閥門的可靠性,對航天閥門的溫度、振動和動作試驗數據進行預處理,融合成統一的先驗分布,參與閥門多源信息融合的可靠性評估過程。該方法利用Bayes理論,對現場數據的評估結果進行修正,減少對試驗樣本量的需求。
機械中的疲勞強度、疲勞壽命、磨損壽命、腐蝕壽命大多服從Weibull分布。根據工程經驗,閥門的可靠工作次數為壽命型參數,其可靠性分布規律符合Weibull分布。因此,本文假設閥門服從兩參數的Weibull分布——(,),其中為形狀參數,為尺度參數(也稱特征壽命)。航天閥門屬于高可靠性產品,受到研制時間和研制經費的限制,研制試驗常常采用截尾試驗,試驗終止時僅出現極少失效甚至無失效。
閥門的溫度試驗采用不等定時截尾的方式,在試驗進行到截尾時間t時,進行第次截尾,去掉t個樣本。已知投入試驗的8個樣本在截尾時均未失效,截尾時間與對應的無失效樣本量t如表1所示。

表1 溫度試驗的不等定時截尾數據Tab.1 Random censored data from temperature test
一組不等定時截尾數據可以看作多組定時截尾數據。采用傳統統計方法中的極大似然估計(MLE)法進行可靠性評估時,若樣本數目足夠多,收斂性質會表現良好,且此方法簡單、結果無偏或漸近無偏。但是,MLE法并不適用于無失效的情況。王玲玲等提出修正MLE(MMLE)法,引入修正常數(0<<1)作為單組無失效數據的偽失效數,再采用MLE方法來開展無失效數據的可靠性評估方法。因此,溫度試驗數據的似然函數為

(1)
式中:和為溫度試驗數據的Weibull參數;為定時截尾組數。
形狀參數和特征壽命的MMLE法估計結果如(2)式所示:

(2)



(3)
顯著性水平可以依據可靠性要求取05、04或03,則特征壽命的評估結果為

(4)
根據相似產品或歷史數據,結合專家經驗給出形狀參數的取值范圍,通過數值迭代法確定(2)式中和
振動試驗以服從高斯分布的隨機振動形式進行,樣本在不同的功率譜密度下,分別進行定時截尾試驗。將6個受試樣本分3組進行恒應力加速壽命試驗,試驗過程中無失效產生,試驗結果如表2所示。已知工況條件下,功率譜密度為2.5/Hz.
在隨機振動下,閥門所受應力與壽命之間的關系滿足逆冪律模型:
=-,
(5)
式中:為振動應力下的閥門壽命;為系統黏性阻尼常數;為振動應力水平;為材料結構常數。

表2 振動試驗加速應力壽命試驗數據Tab.2 Accelerated stress life test data from vibration test
(5)式兩邊取對數,得
ln=0+1ln,
(6)
式中:0和1分別為加速方程系數、常數。
定義環境C對環境B(工況環境)的環境因子為

(7)
式中:為產品在環境B即工況環境下的壽命;為產品在環境C下的壽命。
由(6)式聯合(7)式推導得到加速方程為逆冪律模型下的環境因子:
=exp (1(ln-ln)),
(8)
式中:、分別為產品在環境B和環境C下的振動應力水平。
根據(8)式將表2中的恒應力加速試驗數據換算為工況應力水平下的數據,進一步地,參考溫度試驗數據的MMLE法評估完成振動試驗數據的可靠性評估。
動作試驗同樣采用不等定時截尾的方式,在進行a次開閉后,第次去掉a個樣本。已知投入試驗的10個樣本在截尾時均無失效,截尾次數a與無失效樣本量a如表3所示。

表3 動作試驗的不等定時截尾數據Tab.3 Random censored data from action test
閥門正常開閉次數是衡量閥門可靠性的重要指標。本文將開閉次數處理為連續型數據,以作為現場數據參與多源信息融合的可靠性評估過程。
參考文獻[13]的Bayes可靠性評估方法,首先給出單獨對動作試驗數據進行可靠性評估的結果。不等定時截尾數據同樣可以看作多組定時截尾數據。假設:形狀參數服從均勻分布(,),、分別為形狀參數的上、下限;任務可靠度服從均勻分布(,1),為任務可靠度的下限值。任務工作次數為,則動作試驗數據的似然函數為

(9)
式中:和為動作試驗數據的Weibull參數;為定時截尾試驗的組數;(a)為動作試驗數據下的可靠度函數。
根據表3中的數據,定時截尾試驗組數=4,在任務工作次數時的可靠度函數為

(10)
變換可得

(11)
(11)式代入似然函數,由Bayes定理得到形狀參數和任務可靠度的聯合驗后分布,并提取核函數:

(12)

馬爾可夫鏈蒙特卡洛(MCMC)抽樣方法的基礎理論是馬爾可夫過程。在MCMC方法中,為在一個指定的分布上抽樣,根據馬爾可夫過程,首先從任一狀態出發,模擬馬爾可夫過程,不斷進行狀態轉移,最終收斂到平穩分布。Metropolis-Hastings(M-H)抽樣算法和Gibbs抽樣算法是MCMC方法中使用較為廣泛的兩種形式。M-H抽樣算法不要求已知分布是對稱的,且算法所構造的馬爾可夫過程滿足細致平穩條件。因此考慮使用M-H抽樣算法計算Bayes驗后分布。
航天閥門具備所有航天產品的共性,可投入試驗的樣本量很小、試驗經常出現無失效。充分利用研制試驗數據,可在一定程度上增加可靠性評估的準確性。在融合壽命相關數據的場合,Bayes方法應用最多, 因為Bayes方法結合現場數據與各類先驗信息,可以到達修正參數估計結果的效果,同時需要的樣本量相對MLE方法來說也更少。將不同類型試驗的數據轉化為先驗分布,與現場數據進行一致性檢驗,融合成統一的先驗分布,可以有效提高閥門Bayes可靠性評估的精度。
試驗數據的預處理工作主要是對數據的可靠性評估。將溫度和動作試驗數據的預處理結果,即可靠度,作為先驗信息,可進一步轉化為先驗分布。
已知失效率的共軛先驗分布為Gamma分布,根據變換下的不變性原則可得,閥門任務工作時間(min)時的可靠度先驗分布為負對數Gamma分布。取任務可靠度的先驗分布為

(13)
先驗分布(13)式中存在兩個超參數和,且>0,>0本文采用最大熵法利用可靠度點估計值確定多源信息的先驗分布參數。
若已知在任務時間時的可靠度點估計值為,使用最大熵法(14)式可以確定負對數Gamma分布(13)式中的超參值:

(14)
進一步簡化(14)式,可得

(15)

這樣,問題就被轉化為一維規劃問題,有利于未知參數、的求解。
使用Bayes方法進行信息融合,首先要進行一致性檢驗,確保所有先驗信息與現場數據是屬于同一整體的。 秩和檢驗、Mood檢驗和假設檢驗都可以直接利用樣本數據進行檢驗,但是要求有一定的失效數,這些方法對無失效數據是行不通的。因此,考慮Bayes置信區間法進行檢驗,由單個信息源結合現場數據得到參數的置信區間,如果現場數據在無信息先驗下的參數點估計值落在置信區間內,則認為該信息源與現場數據是相容的,即通過一致性檢驗。


(16)


為融合多源信息,需要先將可靠度的先驗分布(13)式轉化為失效率的先驗分布,可靠度函數為
=exp (-)
(17)
對可靠度函數關于失效率求導,可得

(18)
(18)式中,導函數恒小于0,因此可靠度函數是關于的減函數,有

(19)
等式兩邊關于求導,可得

(20)
假設形狀參數服從均勻分布(,),其中,參數上、下限和根據相似產品、歷史數據由專家經驗給出,則形狀參數與失效率的聯合先驗分布為

(21)
21節解出的超參數和代入(21)式,分別得到溫度和振動試驗信息所形成的先驗分布(,)和(,),則融合后的先驗分布為
(,)=(,)+(,),
(22)
式中:、分別為溫度試驗、振動試驗的先驗信息在融合中所占的權重。根據第Ⅱ類極大似然估計(ML-Ⅱ)法可以確定權重值。
假設某現場數據樣本服從分布(|),為未知參數。已知多源試驗信息形成的先驗分布為ρ(),其中=1,2,…,由先驗分布ρ()得到的邊緣分布為

(23)
將(|)看作是的函數(極大似然函數),它與的先驗分布無關,只反映了樣本的信息。ML-Ⅱ法實質上就是根據在不同的先驗分布下,現場樣本出現的似然性大小來確定不同的融合權重。
假設現場數據,,…,,將其看作是由邊緣分布(,ρ)產生的,得到類試驗先驗分布下的現場數據似然函數為

(24)
則權重的計算方法為

(25)
動作試驗數據作為現場數據,則閥門的極大似然函數為
(|,)=∏(a;,)=exp (-),
(26)

依據貝葉斯原理,融合后的驗后分布為

(27)
驗后分布計算使用MCMC抽樣方法中的M-H算法,采用MATLAB軟件實現的抽樣過程如圖1所示。
已得形狀參數與失效率的點估計值,由(28)式計算得到Weibull分布中特征壽命的估計值,至此完成航天閥門多源信息融合的可靠性評估。

(28)
溫度試驗數據依據11節方法,給定形狀參數的區間范圍為[1,6],數值迭代法得到溫度試驗下的閥門Weibull分布參數:=327,=70724 6

根據表2中的加速應力(振動)試驗數據擬合加速模型,獲得加速模型系數:
ln=4470 3-0484 8ln,
(29)
使用環境因子將加速應力數據折算到工況水平,結果如表4所示。

表4 振動試驗數據折算結果Tab.4 Conversion results of vibration test data
折算后的數據仍具備無失效特性,MMLE方法進行參數估計得到振動試驗下的閥門Weibull分布參數:形狀參數=3,特征壽命=154.662 1.
單獨對動作試驗無失效數據進行可靠性評估時,采用Bayes方法,給定形狀參數取值范圍為 [2,3],任務工作次數為20次,任務可靠度下限為0.8,得到動作試驗下的閥門Weibull分布參數:=2.284 7,=831.818 9.
假設閥門的任務工作時間為10 min,動作試驗任務開閉次數為20次。給定的形狀參數取值范圍為[2,3]。根據以上可靠性評估結果可以得到各類試驗的先驗分布參數,如表5所示。
3類試驗形成的先驗信息已知,對其進行一致 性檢驗,顯著性水平取0.1,計算得到溫度、振動試驗數據下的可靠度區間估計和動作試驗作為現場數據的無信息先驗點估計值,如表6所示。

表5 3類試驗的先驗分布參數Tab.5 Prior distribution parameters for three types of tests

表6 一致性檢驗計算結果Tab.6 Calculated results of consistency test
現場數據的點估計值落入區間估計范圍內,滿足Bayes置信區間法對多源信息通過一致性檢驗的要求。因此,3類試驗信息可以進行多源信息融合。結合表1與表4中的數據,依據(24)式和(25)式計算得到融合權重為==0.5.
經推導,驗后分布的核函數為
(|)∝()(|)=
((,)+(,))exp (-)
(30)
通過MCMC抽樣方法得到閥門分布參數的點估計值為=1595 355×10、=255,抽樣迭代結果頻數分布分別如圖2和圖3所示。將抽樣得到的參數估計值代入(28)式,得特征壽命=466007 7

圖2 MCMC抽樣失效率λ頻數分布圖Fig.2 Frequency distribution of failure rate λ by MCMC sampling method

圖3 MCMC抽樣形狀參數m頻數分布圖Fig.3 Frequency distribution of shape parameter mby MCMC sampling method
為了更加直觀地展示可靠性評估的結果,將3類試驗數據單組的可靠性評估結果與融合后的可靠性評估結果歸納于表7,并給出任務可靠度,在圖4中畫出對應的可靠度函數曲線。

表7 不同類型試驗數據和多源信息融合的 可靠性評估結果Tab.7 Reliability evaluation results of different types of test data and multi-source information fusion

圖4 可靠度函數曲線Fig.4 Reliability curves
已知某型號閥門進行3個批次的生產,共組裝37個閥門。這些閥門經過了典型試驗,包括常溫性能檢查試驗、高低溫性能檢查試驗、振動后性能檢查試驗和動作壽命試驗,又進行了延壽試驗。在系列試驗中共出現3個失效,失效時的累積動作次數分別為2 050、5 450和4 400. 為了簡化計算過程,假設該批次閥門的壽命服從指數分布。通過MLE法計算任務動作次數為20次時的可靠度為0.995.
本文方法計算出的任務可靠度0.999 7與大樣本極大似然的計算結果0.995相比較,誤差在0.5%以內。此外,在第2 050次動作試驗失效的閥門實際上經歷了1 000次的高溫動作壽命試驗,因此,閥門的真實任務可靠度大于0.995. 閥門設計要求為0.999 83,專家一致認為該型號閥門達到任務可靠度指標。通過以上對比可以有效說明本文評估方法的準確性和精確度。
為了進一步驗證本文方法的穩定性,計算機仿真運行10 000次,計算任務可靠度的平均值和平均相對誤差,結果如表8所示。

表8 方法穩定性驗證結果Tab.8 Stability verification results of the proposed method
由表8可知,平均相對誤差相較可靠度平均值而言非常小,說明本文方法具有較好的穩定性。
在航天閥門研制過程中,定量的可靠性評估還缺乏理論方法的支持。航天閥門投入可靠性試驗的樣本量小,而且屬于高可靠性產品,在試驗過程中往往僅能得到有限的無失效試驗數據。
針對航天閥門不同試驗數據的特點,需要采用合適的方法對數據分別進行預處理。將溫度試驗和振動試驗數據信息融合為統一的先驗分布,動作數據作為現場數據,基于Bayes原理完成對航天閥門的可靠性評估。基于多源信息融合的航天閥門可靠性評估方法,既解決了小子樣、無失效數據的可靠性評估問題,又充分利用已有數據,提高了閥門可靠性評估的精度和可信度。