基于混合Gamma分布的機載產(chǎn)品可靠性建模

李軍亮, 祝華遠, 王 正, 王利明, 張鑫磊

(海軍航空大學(xué)青島校區(qū), 山東 青島 266041)

0 引 言

艦載機在其壽命周期內(nèi),要遂行多種作戰(zhàn)任務(wù),需在沿海、近海、遠海等環(huán)境中服役,機載設(shè)備受到溫度、濕度、氣壓、鹽霧、振動、加速度等多種自然環(huán)境應(yīng)力和平臺環(huán)境應(yīng)力耦合作用,故障率高且失效模式復(fù)雜,嚴重影響艦載機作戰(zhàn)效能的發(fā)揮。因此建立科學(xué)的評價方法,并準確評價艦載機機載產(chǎn)品的可靠性是亟待解決的問題。

機載產(chǎn)品主要包括機械類、電子類以及機電耦合類產(chǎn)品,具有小批量、多批次、高可靠、長壽命等特點,研制生產(chǎn)階段的可靠性試驗只能收集在特定應(yīng)力下的產(chǎn)品可靠性信息,不能完全模擬真實使用環(huán)境,且GJB-899A主要對壽命分布為指數(shù)分布的設(shè)備級產(chǎn)品的可靠性試驗進行規(guī)范,未能考慮非指數(shù)型分布或具有性能退化特征的產(chǎn)品特征,適用范圍有限[1]。工程中處理壽命分布的方法有3種：一是根據(jù)先驗信息(主/客觀信息)選擇適用的分布;二是根據(jù)經(jīng)典分布模型進行壽命數(shù)據(jù)擬合、擬合優(yōu)度檢驗以及參數(shù)估計[2-3];三是采用貝葉斯統(tǒng)計方法綜合上述兩種方法。常用的壽命分布模型有指數(shù)、正態(tài)、對數(shù)正態(tài)、威布爾和Gamma分布等[4],但是不同的分布模型都有各自的使用范圍,工程經(jīng)驗表明，電子產(chǎn)品的失效多為指數(shù)分布[5],機械部件的疲勞過程則服從威布爾分布或?qū)?shù)正態(tài)分布[6-7]。機載產(chǎn)品結(jié)構(gòu)復(fù)雜,一般是機械、電子和液壓等器件的組合,在實際使用過程中故障行為復(fù)雜,往往存在故障傳播、故障耦合以及共因失效等特征[8-9],因此單一分布的擬合能力有限。Gamma分布既可以擬合產(chǎn)品的退化數(shù)據(jù),也可以擬合其壽命數(shù)據(jù),黃卓等[10]證明了有限Gamma分布的稠密性,并采用混合Gamma分布構(gòu)建了一種通用的壽命數(shù)據(jù)擬合方法[11],基于最大期望(expectation maximun， EM)算法求解分布參數(shù),并設(shè)計了混合分布的分支控制策略,其擬合結(jié)果與單一分布擬合相比更為精確。在此,擬根據(jù)Gamma分布的自身特性,構(gòu)建一種更為簡潔和有效的混合分布模型和優(yōu)化求解算法,并結(jié)合實際算例對構(gòu)建方法進行驗證。

1 三維混合Gamma 分布的壽命分布擬合

有限混合 Gamma 分布在[0, +∞]上的全體概率分布函數(shù)中稠密[10],即可以逼近任意壽命分布函數(shù)。假設(shè)任意系統(tǒng)故障時間分布函數(shù)可由多個Gamma分布構(gòu)成,不同產(chǎn)品的故障密度概率可表示為

(1)

分析式(1)可知,要想明確系統(tǒng)的壽命分布曲線,首先要確定式(1)中的各項參數(shù)集{al,αl,λl}以及M,文獻[11]在其推理過程中并未明確給出M的確定方法,在算例分析中假設(shè)M≤5,和分布參數(shù)一起求解,求解過程比較復(fù)雜。分析Gamma分布的特點可知,當形狀參數(shù)αl∈(0,+∞)時,其故障率特點可以分為3種趨勢:單調(diào)增加、單調(diào)遞減和恒定不變。當形狀參數(shù)0<α<1,α=1和α>1時(λ=1),分別代表故障率遞減、恒定或者遞增,如圖1所示。

圖1 不同參數(shù)的Gamma分布形狀 Fig1 Shape of Gamma distribution with different parameters

又因為有限維混合Gamma分布在實數(shù)概率空間中稠密[10-11],理論上不同形狀的分布參數(shù)可以通過智能優(yōu)化算法求解,從而獲得最優(yōu)解。因此,擬對M的混合策略進行簡化,假設(shè)任意壽命分布曲線可由3條Gamma曲線組成,即maxM≤3,然后根據(jù)優(yōu)化模型的優(yōu)化目標精度控制擬合精度,從而控制各支分布的參數(shù),其概率密度函數(shù)為

(2)

其故障分布函數(shù)為

(3)

定理1由3支Gamma 分布組成的混合Gamma 分布在[0, +∞) 上的全體概率分布函數(shù)中稠密。

證明假設(shè)g(x)為任意的概率密度函數(shù),其數(shù)學(xué)期望為E(X)=θ,其方差為D(X)=σ2,欲證明f(x)在g(x)中稠密,只需證明存在度量函數(shù)ρ和ε→0,使得ρ≤ε,則原假設(shè)成立。

兩邊取平方可得

證畢

根據(jù)定理1可知,由3支Gamma分布組成的混合分布在[0, +∞)的故障分布函數(shù)中稠密,即式(3)可以擬合出任意故障分布函數(shù),且參數(shù)集中的未知參數(shù)集為{al,αl,λl},而文獻[11]的算例中則含有15個未知參數(shù),因此本文的混合模型未知參數(shù)更少,求解過程相對簡單。特別的,當al=1時,則混合分布為Gamma分布[13-14]。

2 三維混合Gamma分布優(yōu)化模型構(gòu)建

定理1雖然證明了混合Gamma分布的稠密性,但是在實際應(yīng)用中分布參數(shù)集都是未知的,只能通過產(chǎn)品在服役環(huán)境下產(chǎn)生的故障樣本得出其樣本分布的一些統(tǒng)計特性的

2.1 小樣本數(shù)據(jù)的優(yōu)化模型

當樣本量n≤20時,一般采用海森公式、數(shù)學(xué)期望公式或者近似中位秩公式來直接代替經(jīng)驗分布函數(shù)值,在此選用中位秩公式來替代經(jīng)驗分布,如下所示:

(4)

式中:n為樣本容量;i為樣本失效時間由小到大排列的順序,即t1≤t2≤…≤tn。

采用均方誤差MSE作為度量測度,即

(5)

當樣本已知時,式(5)中含有未知參數(shù)集{al,αl,λl},優(yōu)化模型如下:

(6)

(7)

當樣本已知時,可以確定分布參數(shù)置信上下限。

2.2 大樣本數(shù)據(jù)的優(yōu)化模型

Person檢驗是一種比較通用的擬合優(yōu)度檢驗方法,Person統(tǒng)計量適用于一般分布,而不針對特定分布。當樣本量較大時,將Person統(tǒng)計量作為擬合模型和真實分布量的測度。

外場數(shù)據(jù)采集的時效性,不可能像在實驗室環(huán)境下一樣精確,尤其是整個機群中同類產(chǎn)品數(shù)據(jù)的采集可能會存在時間誤差和隨機誤差,為了減少這種誤差的影響,可以根據(jù)樣本情況對采集的時間區(qū)間進行劃分。在此,根據(jù)樣本量n將樣本按時間取值范圍分為k個時間區(qū)間,則ni為第i個區(qū)間內(nèi)的失效樣本數(shù),總體X落入第i個區(qū)間的概率pi為

pi=F(t(i))-F(t(i-1)),i=1,2,…,k

(8)

則Person檢驗統(tǒng)計量可表示為

(9)

將式(3)代入式(9)可得

χ2=

(10)

構(gòu)建以下優(yōu)化模型:

minχ2=

(11)

(12)

3 基于改進粒子群算法的優(yōu)化算法設(shè)計

單一分布模型中比較常見的分布參數(shù)評估方法有最小二乘法、極大似然估計(maximum likelihood estimate, MLE)、EM算法等。文獻[11]在EM算法的基礎(chǔ)上設(shè)計了控制分支M的一種迭代算法,用來求解有限混合Gamma分布的參數(shù)求解,其求解方法和過程都比較復(fù)雜。近年來,遺傳算法、粒子群算法等群智能算法被用來求解可靠性參數(shù)優(yōu)化問題,并且取得了不錯的效果[14-15]。其中,粒子群算法收斂速度快,且模型構(gòu)建過程比較簡單,在此采用自適應(yīng)粒子群算法對混合分布的參數(shù)進行求解,算法的基本流程設(shè)計如圖2所示。

圖2 基于自適應(yīng)粒子群優(yōu)化算法的參數(shù)求解流程Fig.2 Parameter solving process based on adaptive particle swarm optimization algorithm

如圖2所示,在使用粒子群算法的時候,首先要確定參數(shù)集P(N,c1,c2,xLIMT,vLIMT,wLIMT,M,D),其中N為群粒子數(shù)目,c1和c2為加速系數(shù),D為種群維度即目標函數(shù)中自變量個數(shù),M為迭代次數(shù),xLIMT和vLIMT分別為根據(jù)約束條件設(shè)定的粒子自變量取值范圍和粒子迭代移動幅度參數(shù),vLIMT為加權(quán)系數(shù)閾值。則在第t時刻,各個粒子位置參數(shù)為Xi(t)={xi,1(t),xi,2(t),…,xi,n(t)},速度公式為Vi(t)={vi,1(t),vi,2(t),…,vi,n(t)},個體最優(yōu)解為pBesti(t)={pBesti,1(t),pBesti,2(t),…,pBesti,n(t)}，全局最優(yōu)解為gBesti(t)={gBesti,1(t),gBesti,2(t),…,gBesti,n(t)},位置和速度的更新方程公式分別為

vi,j(t+1)=wvi,j(t)+c1r1[pi,j-xi,j(t)]+c2r2[pg,j-xi,j(t)]

(13)

xi, j(t+1)=xi, j(t)+vi, j(t+1),j=1,2,…,d

(14)

式中:r1和r2分別為隨機數(shù)介于0～1。加權(quán)系數(shù)的大小會影響算法的收斂性能,自適應(yīng)權(quán)重可以有效避免陷入局部最優(yōu)解,從而求得全局最優(yōu)解,自適應(yīng)權(quán)重更新公式為

(15)

式中:fmin和favg分別為最小和平均目標值；wmin和wmax為最小和平均權(quán)重。

4 算例分析

不同類型的機載產(chǎn)品的失效機理不同,在不同環(huán)境表現(xiàn)出的故障行為也不一致。某艦載直升機機群首次大修間隔內(nèi)故障統(tǒng)計數(shù)據(jù)表明,機載電子產(chǎn)品和機電產(chǎn)品的故障率相對較高,本節(jié)分別以某型艦載直升機多功能顯示器和微動開關(guān)的故障樣本為例來驗證算法的正確性。

4.1 算例1

多功能顯示器故障主要是因為其內(nèi)部的印制電路板(printed circuit board, PCB)失效引起,在其服役過程中受到隨機振動、酸性氣體、鹽霧以及溫濕度的作用,該設(shè)備的密封結(jié)構(gòu)、防腐涂層失效進一步引起PCB上的叉指電極、通孔等發(fā)生腐蝕等故障行為。其樣本故障如表1所示,表中第5列Fn(ti)為采用近似中位秩公式計算的壽命分布近似值。采用本文設(shè)計方法,將計算結(jié)果Fn(ti)代入式(6)并展開可得優(yōu)化目標函數(shù)。

表1 某設(shè)備PCB板故障數(shù)據(jù)樣本

(16)

分別采用經(jīng)典粒子群算法和自適應(yīng)粒子群算法進行求解,其中粒子群算法參數(shù)為算法的初始參數(shù)P(50,3,3,0.8,100,9),自適應(yīng)粒子群算法的參數(shù)設(shè)置為P(50,3,3,0.8,0.6,100,9),其計算結(jié)果如圖3所示。如圖3所示,兩種算法收斂速度較快,且適應(yīng)度函數(shù)均接近0.003 771 096 562 428,其中自適應(yīng)用粒子群算法的參數(shù)集的優(yōu)化結(jié)果為(0.9, 12.19, 221.12, 0.09, 1, 221.12, 0.01, 0.2, 221.12),將該參數(shù)代入式(3)可得產(chǎn)品的壽命分布函數(shù),如下所示:

F(t)=0.9Γ(t,12.194 1,221.12)+0.09Γ(t,1,1.03)+ 0.01Γ(t,20,0.99)

(17)

圖3 不同算法的收斂性及最小適應(yīng)值Fig.3 Convergence of different algorithms and fitness value

為了更直觀地顯示壽命分布的擬合曲線,分別將混合分布模型與近似中位秩法以及殘存比率法進行比較,結(jié)果如圖4所示,系統(tǒng)的可靠度如圖5所示。如圖4所示,近似中位法計算的故障率最高,殘存比例最低,混合分布的故障函數(shù)值介于近似中位秩和殘存比例法中間,因為近似中位秩法不考慮刪失樣本數(shù),所以數(shù)值偏大,而殘存比率法需要考慮刪失樣本,所以數(shù)值偏小。實際上,近似中位秩法更適用于完全樣本的故障率近似,但是在實際使用過程中很難獲得產(chǎn)品的完全樣本,尤其是航空產(chǎn)品具有高可靠、長壽命的特征,在一定時間內(nèi)很難獲得產(chǎn)品的完全樣本,其樣本一般為截尾數(shù)據(jù),殘存比率法在工程中經(jīng)常被用于計算具有刪失特性的截尾樣本數(shù)據(jù),但其計算結(jié)果偏于保守,因此混合Gamma分布擬合方法可以有效改良兩種方法的缺陷,接近產(chǎn)品的真實故障分布特性。

圖4 3種算法的比較Fig.4 Comparison of three algorithms

圖5 PCB板的可靠度函數(shù)Fig.5 Reliability function of PCB board

在裝備保障的實際過程中,可結(jié)合式(17)和圖5合理設(shè)定產(chǎn)品的可靠性閾值或者剩余壽命,從而有針對性地制定預(yù)防性維修策略。

4.2 算例2

以機載微動開關(guān)為研究對象,在艦載服役環(huán)境中,微動開關(guān)易受到酸性鹽霧及溫濕度的影響,開關(guān)內(nèi)部觸點發(fā)生腐蝕從而引起開關(guān)接觸電阻變大,直至功能失效,因此其失效過程為退化失效和突發(fā)失效的綜合作用。其失效前后的照片如圖6所示,圖6(a)為裝機前的新樣品,圖6(b)為失效后拆卸下來的樣本。

圖6 DK1-2微動開關(guān)的失效前后照片F(xiàn)ig.6 Photos of DK1-2 microswitch before and after the failure

假設(shè)機群中每架飛機的使用環(huán)境和使用頻次基本一致,采集其近5年的故障數(shù)據(jù),得到80組故障數(shù)據(jù),如表2所示。

表2 微動開關(guān)故障數(shù)據(jù)樣本

如表2所示,樣本數(shù)據(jù)中最小故障時間為218 h,最大故障時間為266 h,樣本容量為80,取時間間隔Δt=10 h,可得樣本的故障頻率直方圖如圖7所示,即可獲得的ni。觀察圖7可知，樣本數(shù)據(jù)非正態(tài)、指數(shù)分布等單一分布,因此按照式(10)構(gòu)建樣本的優(yōu)化目標函數(shù)為

minχ2=

(18)

圖7 微動開關(guān)故障樣本頻率分布直方圖Fig.7 Histogram of frequency distribution of micro switch fault samples

采用文中設(shè)計方法求解可得,參數(shù)設(shè)置為(500, 3, 3, 0.8, 0.5, 200, 9),算法收斂過程如圖8所示,最佳適應(yīng)值為0.001 057 3,優(yōu)化結(jié)果的參數(shù)集為(0.95, 260, 1, 0.02, 1, 1.03, 0.03, 30, 0.99)。

圖8 算例2的收斂過程Fig.8 Convergence process of Example 2

將所求的參數(shù)代入式(3)可得產(chǎn)品的壽命分布模型為

F(t)=0.95Γ(t,260,1)+0.02Γ(t,1,1.03)+ 0.03Γ(t,30,0.99)

(19)

可靠度曲線如圖9所示。工程實踐中,一般以威布爾分布或者對數(shù)正態(tài)分布擬合機電產(chǎn)品的故障曲線[13-14],圖9中“△”的曲線為以算例2的樣本按照99.5%置信度擬合出的服從對數(shù)正態(tài)分布的產(chǎn)品可靠度曲線;“*”為采用本文設(shè)計的算法計算的產(chǎn)品可靠度;“+”為采用文獻[11]設(shè)計的算法計算的產(chǎn)品的可靠度。其中,“△”曲線在使用時間小于219 h之前,認為產(chǎn)品可靠度一直大于99%,在使用時間接近220 h產(chǎn)品可靠度迅速衰減,接近0,與表2中的故障樣本數(shù)據(jù)反映的特征有所差別,不能有效反映出產(chǎn)品的真實可靠性和故障特性。本文設(shè)計算法和文獻[11]中的算法計算的可靠度曲線形狀比較相似,但文獻[11]中的算法計算產(chǎn)品可靠度略低于本文算法,本文算法的可靠度值基本介于對數(shù)正態(tài)分布和文獻[11]算法之間,且與表2中的故障數(shù)據(jù)比較吻合,可以較為有效表現(xiàn)產(chǎn)品的可靠性特征。

圖9 微動開關(guān)可靠度Fig.9 Product reliability of micro switch

如圖9所示,微動開關(guān)的可靠性在不同時間內(nèi)呈現(xiàn)出不同的特征,在0～40 h產(chǎn)品可靠性較高,但處于下降階段;在40～227 h之間可靠度趨于平穩(wěn),其最小值大于0.94;在270 h之后可靠度迅速下降;在300 h之后產(chǎn)品的可靠度趨近于零。因此,可結(jié)合式(19)在不同的壽命階段合理設(shè)定產(chǎn)品的可靠性閾值以及預(yù)防性維修策略,將比目前定期檢查策略更為有效[16-18]。

5 結(jié) 論

本文分析了在不需要借助工程經(jīng)驗判斷產(chǎn)品壽命分布模型類型的前提下,如何基于故障數(shù)據(jù)構(gòu)建產(chǎn)品的混合Gamma分布模型,并從理論上證明了該方法的可行性,設(shè)計了求解模型的粒子群算法,并通過算例進行了驗證和分析。算例1結(jié)果顯示,當樣本量小于20時,本文方法與樣本真實誤差小于0.000 3,且分別與近似中位秩法和殘存比例法進行了比較,證明了本文方法更為精確。算例2說明在大樣本情況下,本文方法同樣適用,最佳適應(yīng)值為0.001 057 3,通過算法和誤差控制可以獲得分布參數(shù)精確數(shù)值解,擬合出合理的產(chǎn)品壽命分布模型。

本文方法適用于在復(fù)雜環(huán)境中服役的復(fù)雜產(chǎn)品,如機電組合設(shè)備,因為其失效機理較為復(fù)雜,受到多種因素的綜合作用,且電子產(chǎn)品和機械設(shè)備的壽命分布模型往往不一致,單一分布模型在擬合其壽命分布模型的精度有限,而本文方法通過混合分布可以更加真實地反映產(chǎn)品的故障分布統(tǒng)計特征。

文中只討論了基于故障數(shù)據(jù)即壽命數(shù)據(jù)的可靠性評估方法,對于系統(tǒng)的退化數(shù)據(jù)并未進行討論,而基于退化數(shù)據(jù)和壽命數(shù)據(jù)相融合的方法更能準確的評估系統(tǒng)的可靠性。另外,艦載產(chǎn)品的服役環(huán)境非常復(fù)雜,目前尚未有有效的艦載機機載產(chǎn)品在服役過程中適用于環(huán)境因子折算的環(huán)境剖面劃分方法以及混合分布模型的環(huán)境因子折算方法,開展此類研究將非常有意義。