基于貝葉斯的飛機部件平均故障間隔的估計

余 芬，張 鑫

(中國民航大學 航空工程學院，天津 300300)

0 引言

改裝是指在航空器及其部件交付后進行的超出其原設計狀態的任何改變，不同方式的飛機改裝中，服務通告(Service Bulletin，SB)改裝的方案優勢更大，同時其帶來的改裝成本也相較于其他改裝途徑要高。對于非強制性的服務通告類型，航空公司需根據自己的實際情況決定是否依服務通告的要求對飛機實施相應的改裝活動，即進行改裝評估。飛機改裝效果的好壞，取決于可靠性指標的估計。服務通告中的一個基本參數就是平均非例行拆換間隔(Mean Time Between Unscheduled Removals，MTBUR)，該參數估算的準確性不僅影響可靠性和維修性，更重要的是關系到服務通告經濟性決策的正確性。但由于缺乏相應的運行和使用數據，使得對參數的評估往往僅由專家經驗得到，如何在數據相對稀缺的環境下獲得準確的估計值成為了服務通告評估的重要內容。

目前，國內相關學者對服務通告的研究大都在服務通告管理策略與控制方法等方面[1-4]，對定量化改裝評估的研究較少。文獻[5]采用極大似然方法估算平均故障間隔時間(Mean Time Between Failure，MTBF)，分析執行服務通告后部件壽命變化對航空公司經濟性的影響。文獻[6]利用凈現值法對服務通告執行的經濟性進行量化。文獻[7]針對壽命服從指數分布的部件應用貝葉斯方法進行可靠性指標參數估計。上述文獻雖對小樣本環境下的飛機改裝提出了一些方法，但大多數飛機部件的壽命并不服從指數分布，對于壽命分布服從非指數分布類型的部件未做研究，同時未考慮改裝樣本量的經濟性問題。

本文以服從雙參數威布爾壽命分布的飛機結構件為研究對象，采用貝葉斯方法對可靠性指標進行參數估計[8]。根據馬氏鏈蒙特卡洛方法，建立后驗平穩分布的馬爾科夫鏈[9-10]，對分布參數的貝葉斯估計進行求解，以解決數值積分問題，保證可靠性評估的有效實施。

1 可靠性指標的貝葉斯估計

1.1 貝葉斯估計基本理論

貝葉斯方法有效地利用了總體、樣本和先驗信息，貝葉斯基本公式為：

(1)

其中：X為采樣樣本；θ為總體分布的參數；Θ為參數空間；π(θ)為參數θ的先驗密度函數；π(θ|X)為參數θ的驗后分布密度函數；L(X|θ)為似然函數。

1.2 參數估計模型的建立

MTBUR是服務通告經濟性評估的重要參數，它被用來計算每飛行小時的成本差異和評估機隊維修水平。MTBF是指系統或部件相鄰兩次故障的平均時間間隔，它是相鄰兩次非例行拆換之間的平均工作時間。一般來說MTBUR比MTBF要小。為準確地對上述參數進行估計現做如下假設：①所有的拆換都是確定失效后的部件更換；②部件壽命分布服從兩參數的威布爾分布；③改裝件的連續運行時間T服從威布爾分布。則其壽命分布密度函數f可表示為：

(2)

其中：λ為比例參數；k為形狀參數；t為部件連續運行時間?，F把MTBUR表示為MR，改裝樣本量充足的情況下，部件的壽命分布函數可表示為:

(3)

其中：m為形狀參數。則樣本似然函數L可以表示為：

(4)

其中：n為部件個數；t1，t2,…，tn為部件壽命值。此時，令參數m、MR的先驗分布分別為伽馬分布和逆伽馬分布，則參數的聯合先驗分布為:

(5)

其中：a、u為形狀參數；b、v為尺度參數；Γ(a)、Γ(u)為伽馬函數在a、u處的取值。通過貝葉斯估計公式進而得到后驗分布函數π：

(6)

其中：T=(t1，…,tn)為總體部件壽命。由于后驗分布的形式復雜，此處借助于馬爾可夫鏈蒙特卡洛方法解決后驗求解問題。

1.3 馬爾可夫鏈蒙特卡洛模擬

在貝葉斯參數估計中，采用蒙特卡洛抽樣方法計算感興趣的量。令p(x,y)=q(x,y)α(x,y)，表示從建議分布q(*|x)產生一個潛在的轉移x→y，然后依概率α(x,y)決定是否轉移，而p(x,y)為馬氏鏈一個轉移核。選取均勻分布為建議分布，則概率α(x,y)取為：

(7)

其中：my和MRy為一次迭代中從建議分布中抽取的建議值；mx和MRx為迭代的初始值。

2 實例仿真分析

假定適用于 B737-300機型的某份服務通告給出，某飛機結構件的MR大約為490 FH～510 FH(飛行小時)之間?，F有一家航空公司擁有100架該機型的機隊，通過對廠家下發的服務通告的分析，擬選取機隊中的一部分飛機實施試裝?，F工程技術部門給出MR的參考值為500 FH，那么如何選取改裝件數做改裝試驗才能準確評估該參數。假設參數MR=500 FH，m=1.8,m、MR的先驗分布分別服從伽馬和逆伽馬分布。則其分布函數可以表示為：

(8)

(9)

(10)

(11)

根據式(6)～式(11)獲得參數的先驗分布，然后根據貝葉斯理論，通過馬爾可夫鏈蒙特卡洛模擬求參數的后驗分布，m和MR的后驗密度分別為π(m|MR，T)、π(MR|m,T)，其分布圖如圖1和圖2所示，m和MR的全部迭代軌跡如圖3和圖4所示。

圖1參數m的后驗密度分布圖圖2參數MR的后驗密度分布圖圖3m的迭代軌跡

圖4 MR的迭代軌跡

由圖1～圖4可以看出，馬爾科夫鏈達到了穩定狀態，參數的平穩分布即為其后驗分布，說明本方法在做威布爾分布參數估計時是有效的。利用MATLAB仿真產生不同改裝樣本量下的參數m、MR的貝葉斯估計見表1、表2。

表1 不同改裝樣本下m的點估計

表2 不同改裝樣本下MR的點估計

由表1、表2數據可以看出，不同樣本數量下的參數估計值誤差不同。數據稀缺環境下的參數估計，隨著樣本數量的增加，后驗期望的波動幅度隨之減小，參數估計更趨近于理想值，更符合實際情況。同時從估計誤差的角度考慮樣本數量變化的影響，顯然單次實驗的誤差是隨機的，但是誤差發展的趨勢是必然的。估計誤差隨著樣本量的增加逐漸減小，我們可以設定一個誤差范圍來根據誤差的發展趨勢反推我們所需的樣本量。在給定估計誤差的范圍下，可以盡可能地選取少量的改裝樣本件，從而節省經濟成本。與同等條件采用極大似然估計獲得的點估計(見表3)相比，貝葉斯參數估計更加準確。

從表2、表3可以看出，當采用極大似然估計時，在抽取的樣本較少時參數估計的誤差較大。隨著樣本量的增加，兩種方法都能很好地進行估計，差別較小，但是由于融合后驗的貝葉斯估計包含了參數的可靠度的先驗信息，得到的結果更為準確，對于大部件壽命的參數估計更為可靠。

表3 極大似然估計下MR的點估計

3 結語

本文從可靠性角度出發，將評價改裝效果的可靠性指標轉化為壽命分布中的參數，運用貝葉斯方法進行估計。利用馬爾可夫鏈蒙特卡洛模擬方法有效地解決了后驗分布積分復雜帶來的不便。同時通過仿真給出了不同樣本量下的參數估計誤差，其估計結果優于傳統的極大似然估計。

通過設定誤差范圍，可以有效地設定可靠性實驗的樣本數。這對于航空公司在決策飛機改裝試驗時具有工程實際意義。同時，改變MR的大小，非例行工作成本也會跟隨MR做相應的變化。因此客觀準確地對MR做出估算對于航空公司工程人員來說意義重大。