基于混合Gamma分布的機載產品可靠性建模

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

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

0 引 言

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

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

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

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

(1)

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

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

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

(2)

其故障分布函數為

(3)

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

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

兩邊取平方可得

證畢

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

2 三維混合Gamma分布優化模型構建

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

2.1 小樣本數據的優化模型

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

(4)

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

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

(5)

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

(6)

(7)

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

2.2 大樣本數據的優化模型

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

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

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

(8)

則Person檢驗統計量可表示為

(9)

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

χ2=

(10)

構建以下優化模型:

minχ2=

(11)

(12)

3 基于改進粒子群算法的優化算法設計

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

圖2 基于自適應粒子群優化算法的參數求解流程Fig.2 Parameter solving process based on adaptive particle swarm optimization algorithm

如圖2所示,在使用粒子群算法的時候,首先要確定參數集P(N,c1,c2,xLIMT,vLIMT,wLIMT,M,D),其中N為群粒子數目,c1和c2為加速系數,D為種群維度即目標函數中自變量個數,M為迭代次數,xLIMT和vLIMT分別為根據約束條件設定的粒子自變量取值范圍和粒子迭代移動幅度參數,vLIMT為加權系數閾值。則在第t時刻,各個粒子位置參數為Xi(t)={xi,1(t),xi,2(t),…,xi,n(t)},速度公式為Vi(t)={vi,1(t),vi,2(t),…,vi,n(t)},個體最優解為pBesti(t)={pBesti,1(t),pBesti,2(t),…,pBesti,n(t)}，全局最優解為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分別為隨機數介于0～1。加權系數的大小會影響算法的收斂性能,自適應權重可以有效避免陷入局部最優解,從而求得全局最優解,自適應權重更新公式為

(15)

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

4 算例分析

不同類型的機載產品的失效機理不同,在不同環境表現出的故障行為也不一致。某艦載直升機機群首次大修間隔內故障統計數據表明,機載電子產品和機電產品的故障率相對較高,本節分別以某型艦載直升機多功能顯示器和微動開關的故障樣本為例來驗證算法的正確性。

4.1 算例1

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

表1 某設備PCB板故障數據樣本

(16)

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

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

(17)

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

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

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

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

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

4.2 算例2

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

圖6 DK1-2微動開關的失效前后照片Fig.6 Photos of DK1-2 microswitch before and after the failure

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

表2 微動開關故障數據樣本

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

minχ2=

(18)

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

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

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

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

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

(19)

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

圖9 微動開關可靠度Fig.9 Product reliability of micro switch

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

5 結 論

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

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

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