某飛機系統故障率敏感性分析研究

何勃，張文瀚，林健

（成都飛機設計研究所，成都 610091）

飛機在不同的服役階段和使用環境下具有不同的故障率特征[1-2]，根據標準化后的歷史外場使用數據，準確預測不同服役階段和使用環境下的系統故障率，識別出影響故障率大小的敏感因素，可為系統故障預防、故障診斷和故障維修提供依據，對飛機完好率和出勤率的提升具有重要意義。傳統的以使用時間為自變量的故障率分析模型只能給出裝備的理論故障率變化趨勢，而綜合考慮已使用年限和使用環境信息的故障率分析模型能為飛機的故障率預測提供更加可信的依據。

對于我國大部分內陸地區，使用環境對飛機系統故障率的影響主要體現在溫度、濕度和氣壓因素上。溫度因素對飛機系統的故障率影響體現在多個方面，如對于電子產品，高溫環境下元器件損壞、參數漂移的概率增加；對于機械產品，高溫環境下運動機構卡滯、摩擦副過熱的概率增加，低溫環境下系統密封性能下降，溫度劇烈變化帶來的配合關系變化，均會導致系統故障率上升。濕度對飛機系統故障率的影響主要體現在高濕環境下機載產品和機體結構更易發生銹蝕故障，線路及電子元器件凝結冷凝水，造成短路故障率上升。氣壓對飛機的故障率同樣存在影響，低氣壓環境下，發動機起動所需功率增大，剎車系統負荷增加。另外，對于機載制氧制氮系統，氣壓對其故障率也有著重要影響。

飛機系統故障率的影響因素較多，影響關系復雜，且不確定，對飛機系統故障率進行敏感性分析，識別故障率各影響因素的重要度，對飛機的外場使用具有重要指導意義。敏感性分析可分為兩大類：局部敏感性分析和全局敏感性分析[3-4]。局部敏感性可定義為輸出變量對輸入變量在固定位置的偏導，該方法簡單直觀，但當模型非線性較強，且輸入變量維度較高時，可能會得到錯誤的分析結果。全局敏感性分析方法通過探索輸入變量的全部分布區域來分析模型的輸出變化，不依賴于輸入的位置信息，其中Sobol法是一種較為常用的全局敏感性分析方法[5-6]。在利用Sobol 方法解決敏感性分析問題時，根據各影響變量確定樣本的故障率大小則是分析的前提。近年來，代理模型的應用可有效解決該問題。常用的代理模型有響應面、支持向量機、神經網絡和Kriging，其中Kriging 是一種較為精確的插值統計模型，在系統的響應預測中已被廣泛應用[7-9]。

本文對飛機系統的故障率敏感性分析問題進行了研究，綜合考慮飛機服役階段和使用環境，以飛機已使用年限和任務地點的自然環境數據為輸入變量，故障率為輸出變量訓練Kriging，可以根據標準化后的歷史外場數據來預測不同使用年限及使用環境條件下的飛機系統故障率。在進行飛機系統故障率敏感性分析時，通過利用Kriging 對隨機矩陣中各樣本點的故障率進行預測，可以更加高效地計算各輸入變量的敏感性指標，為飛機系統故障預防、故障診斷和故障維修提供依據。

1 方差分解

假設模型輸出響應為Y=y(X)，其中X={X1,X2, …,Xn}，為服從概率密度函數f(X)的n維變量。根據高維模型描述，y(X)可分解為[6,10-11]：

其中：

其中：X～i為除Xi外的其余變量。當X1,X2,…,Xn之間相互獨立時，式（1）等號右邊各項正交。因此，對式（1）等號兩邊取方差可得：

其中：

式（3）等號兩邊同除以V(Y)，可得：

其中：

Si定義為主指標，也稱一階效應指標；Sij=Vij/V(Y)，為二階指標，表征變量Xi和Xj間的二階交叉效應。類似地，s階指標可解釋為s個變量對輸出方差的s階交叉效應。文獻[6]中又提出了全效應指標：

根據全方差公式[12]：

則全效應指標為：

2 Kriging 代理模型

2.1 Kriging 模型的基本思想

Kriging 模型可以作為黑箱函數來預測n維輸入變量X的輸出g(X)[13]，其在敏感性分析中的效率已經被證明[14-15]。Kriging 模型可分解為線性回歸部分和隨機過程[16]：

式中：f(X)=[f1(X),f2(X),….,fm(X)]T稱為基函數；β=[β1,β2,…,βm]T為回歸向量矩陣；fT(X)β表示輸出的平均近似；Z(X)表示回歸模型的偏差模擬。假設為靜態高斯隨機過程，樣本Xk和Xj的協方差為：

其中，k,j=1,2, … ,N，σ2為方差，R(·,·)為高斯相關函數：

式中：θl是相關參數；Xkl和Xjl為向量Xk和Xj的第l個分量。

參數β和σ2可通過式（13）、（14）估計：

式中：F是由f(X)在訓練樣本點的評估值組成的向量；g為對應輸出；R為相關矩陣；Rk,j=R(Xk,Xj)。Rk,j中的未知參數θl可以利用最大似然估計法計算：

Kriging 預測的期望值(X) 為g(X)的最優線性無偏估計：

其中，r(X)是相關向量，r(X)=[R(X,X1), …,R(X,XN)]。

2.2 基于訓練樣本點的Kriging 交叉驗證

當利用Kriging 模型來預測樣本輸出時，模型的構造質量對預測精度起著決定性作用。為了得到更加精確的預測結果，在預測前需對Kriging 模型進行精度驗證，最直接的驗證方法就是隨機抽取若干樣本點作為測試點，通過比較測試點的預測值和真實值之間的誤差來檢驗模型精度，但當測試點的真實輸出值無法直接獲取時，該驗證方法則無法實施。因此，本文引入交叉驗證法充分利用已有信息來驗證Kriging 模型的精度。

交叉驗證的主要思想：對于N個已知真實輸出值gT={g(X1),g(X2), …,g(XN)}的訓練樣本XT={X1,X2, …,XN}，記除去第j個樣本Xj外的其他N-1 個樣本組成的集合為XT/Xj={X1,X2, …,Xj-1,Xj+1, …,XN}?；贜-1 個樣本，XT/Xj與對應真實輸出值gT/g(Xj)訓練的Kriging 模型記為（X），模型（X）在第j個訓練樣本點預測的相對誤差rj為：

依次計算所有訓練樣本的預測誤差，則可以得到當前Kriging 模型的預測誤差：

當利用交叉驗證法驗證Kriging 模型滿足精度要求時，即rRPRESS＜Rs，則認為當前的Kriging 模型可用于預測樣本輸出，否則需再次增加訓練樣本數量來進一步提高Kriging 模型的預測精度。

3 敏感性指標的計算

3.1 直接法

敏感性指標最簡單直觀的計算方法為直接法，首先應用蒙特卡洛法根據輸入變量的概率密度函數f(X)生成Nin組樣本，將每個樣本Xi分別固定在Nin個不同點上。再根據蒙特卡洛法對X～i抽取Nout組樣本，并計算，然后根據Xi的Nin個樣本對應的計算，最終可利用式（6）得到輸入變量Xi的主效應指數。

根據上述分析，當采用直接法來計算時，對于n維輸入變量X1,X2, …,Xn,共需n×Nin×Nout次樣本輸出值計算。假設得到總方差還需2k次樣本輸出計算，則直接法計算主效應指標總共需2k+n×Nin×Nout次樣本輸出值計算。采用直接法計算全效應指標與上述過程相同。

3.2 Sobol 法

當n維輸入變量X1,X2,…,Xn相互獨立時，可以利用Sobol 法來計算各敏感性指標。記A和B為2 個相互獨立的隨機抽樣矩陣：

其中k為抽樣次數，根據文獻[17]和[18]，有：

對于總方差，有：

將式（22）和（23）分別代入式（6）和（9）則可得到輸入變量Xi的主效應指標和全效應指標。

4 故障率敏感性分析方法

裝備故障率主要受服役階段及使用環境因素的影響[19]，在飛機外場保障中，系統的故障率預測對飛機的使用維護具有至關重要的意義。飛機系統的故障率P可根據單位日歷時間內系統發生的故障次數及工作時間計算：

式中：n為單位日歷時間內的系統故障數量；h為單位日歷時間內的系統工作時間。

根據裝備或產品使用規律，其故障率會隨服役階段發生變化。例如，對于大部分機械產品，其故障率趨勢一般呈現浴盆曲線特征，主要分為早期失效期、偶然失效期和耗損失效期3 個階段，所以服役階段為影響裝備故障率特征的一個重要因素。對于我國絕大多數內陸地區，環境因素對故障率的影響主要體現在溫度、濕度和氣壓3 個因素上。原始的使用環境數據為按固定時間間隔采集的一系列數據序列，數據規模龐大，無法直接用于系統故障率敏感性分析。為了提高環境數據處理效率，有必要從環境數據序列中抽取統計特征，首先選取平均溫度、平均濕度和平均氣壓作為故障率敏感性分析的變量，同時根據裝備的外場使用經驗，溫差也會引起系統的故障率變化，平均溫差也應作為故障率敏感性分析的變量。因此，飛機系統的故障率P可表示為：

式中：X1為使用年限；X2為平均溫度；X3為平均溫差；X4為平均濕度；X5為平均氣壓。

由于真實的故障率函數g(X)無法獲得，文中采用Kriging 模型來預測故障率g(X)的數值，使用年限、平均溫度、平均溫差、平均濕度和平均氣壓則作為Kriging 模型的輸入。將Kriging 模型引入Sobol 方法，用Kriging 預測Sobol 法中隨機矩陣各樣本點對應的故障率，并以此來計算各輸入變量的敏感性指標。該方法首先建立初始Kriging 模型，并通過交叉驗證來保證初始Kriging 模型具有一定置信度來預測樣本矩陣A、B、和中各樣本的故障率，進而得到各輸入變量的敏感性指標。綜上，為以更高的計算效率得到滿足要求的故障率敏感性指標，發展了以下的計算程序：

1）分別按時間周期T采集飛機系統在不同使用環境下的環境數據，包括環境溫度、濕度和氣壓數據序列，并統計飛機系統的已使用年限及每個周期T內對應的工作時間和故障數量。

2）計算每個周期T內采集的溫度、濕度和氣壓數據序列的統計特征，即平均溫度、平均溫差、平均濕度和平均氣壓，以及對應的故障率。

3）根據第2 節中Kriging 模型參數建立Kriging模型。相關模型選擇為高斯模型，回歸模型選為常數。文中Kriging 模型借助于DACE 工具箱[20]建立。

4）將不同周期T對應的系統已使用年限、平均溫度、平均溫差、平均濕度和平均氣壓及故障率作為訓練樣本，對Kriging 模型進行訓練。

5）根據式（19）計算當前Kriging 模型的預測誤差rRPRESS。判斷rRPRESS＜Rs是否滿足。若不滿足，則前往步驟1）繼續收集訓練樣本；若滿足，則前往步驟6）。

6）利用蒙特卡洛方法根據使用環境及使用年限的統計分布生成k組樣本點，組成隨機矩陣A。再次利用蒙特卡洛方法根據使用環境及使用年限的統計分布生成k組樣本點，組成隨機矩陣B。

7）將隨機矩陣A和B的所有樣本點（矩陣中的每一行）組成的集合記為St，利用Kriging 模型預測St中每個樣本點對應的故障率。

8）根據Kriging 模型預測的St中各樣本點的故障率計算V(Y)。

9）令i=1，將矩陣和中所有樣本組成的集合記為S(i)，利用Kriging 預測S(i)中每個樣本的故障率。

10）根據Kriging 對S(i)中各樣本的故障率預測結果，利用式（6）和（9）計算變量Xi的主效應指標Si及全效應指標。判斷i=n（n為輸入變量個數）是否滿足，若滿足，則結束算法；否則令i=i+1，回到步驟9）。

飛機故障率敏感性分析方法的流程如圖1 所示。雖然Sobol 方法較直接法在計算效率上已有極大提高，但由于目前積累的外場故障率樣本數量較少，無法滿足直接應用Sobol 方法進行敏感性分析的要求，所以本文選取了訓練代理模型作為黑箱函數預測的方法來對Sobol 方法抽取的已知使用年限和環境變量的樣本進行故障率預測。在代理模型中，Kriging 模型因為計算量相對較小，且在大量文獻[7,9]中被證明對非線性復雜問題等具有較好的近似效果，所以本文選取了Kriging 模型來預測Sobol 方法中所需樣本的故障率。

圖1 飛機系統故障率敏感性分析程序Fig.1 Flow chart of sensitivity analysis for aircraft system failure rate

5 應用研究

本節統計了2016—2018 年某飛機子系統外場不同服役階段和使用環境下對應的系統故障率數據。時間周期T取為日歷月，每個時間周期內按6 h 時間間隔采集溫度、濕度和氣壓數據，根據采集的環境數據序列計算每個周期對應的使用環境統計特征，即平均溫度、平均溫差、平均濕度和平均氣壓。同時統計每個時間周期T內系統故障數量及系統工作時間，進而得到每個時間周期T對應的系統故障率P。以每個時間周期T對應的已使用壽命和使用環境統計特征及其故障率作為訓練樣本，本算例中共計收集了144 組樣本，樣本的故障率數據如圖2 所示。

圖2 交叉驗證過程中樣本點的Kriging故障率預測值與真實值比較Fig.2 Comparison between the value predicted by Kriging and the actual value of sample in the cross-validation process

5.1 故障率預測模型建立

根據144 組訓練樣本及對應故障率建立并訓練Kriging：相關模型選擇為高斯模型，回歸模型選為常數。采用交叉驗證策略（143 組建模，1 組驗證）對初始 Kriging 模型進行精度驗證：計算rRPRESS為0.036，可滿足故障率預測要求。圖2 中對交叉驗證過程中各樣本點的Kriging 故障率預測值與真實值進行了比較。

5.2 故障率影響因素敏感性分析

5.2.1 故障率影響參量的分布

影響飛機系統故障率的變量包括使用年限、平均溫度、平均溫差、平均濕度和平均氣壓，各變量的分布范圍參考系統的實際使用工況給出，分布形式選取為均勻分布，分布參數見表1。

表1 故障率影響因素的分布參數Tab.1 Distribution parameters of the factors affecting failure rate

5.2.2 Kriging 模型對隨機抽樣矩陣的預測

利用蒙特卡洛方法，根據表1 中系統故障率影響變量的分布抽取樣本點，組成隨機矩陣A和B（樣本數選為106），基于5.1 節訓練完成的Kriging 模型預測矩陣A和B組成的樣本集St中每個樣本點的故障率，根據Kriging 對樣本集St的每個樣本的故障率預測值計算總方差V(Y)。完成總方差V(Y)的計算后，繼續利用當前 Kriging 模型預測和(i=1,2,…,5)組成的樣本集S(i)中每個樣本的故障率，然后根據樣本群S(i)中每個樣本的故障率預測值計算及，進而得到該系統故障率5 個影響變量的主效應指標和全效應指標，見表2。

表2 故障率影響變量的主效應指標及全效應指標Tab.2 Main effect and total effect indexes of the failure rate variables

5.2.3 分析結果討論

從表2 中該系統故障率的敏感性分析結果來看，各影響變量的主效應指標與全效應指標的排序較為一致，根據各輸入變量的重要度排序可以定性地比較各輸入變量對系統故障率的影響大小。從主效應指標及全效應指標的排序結果來看，平均溫度和平均濕度為對該系統故障率影響較大的變量。從設計角度來講，要控制系統的故障率，這些變量應當是需要重點關注的。其他的一些變量，如平均溫差和平均氣壓對故障率的影響較小，所以當需要降低故障率計算模型的輸入變量維度或簡化分析模型時，這些變量應是首先可以忽略的。

由于該算例中Kriging 訓練樣本的使用年限分布較多集中于系統的穩定運行期內，系統使用初期及接近壽命極限階段內的樣本量較少，導致Kriging 模型對這2 個階段的信息學習不夠。因此，從本文的敏感性分析結果來看，故障率對使用年限因素不夠敏感，而這可能與實際使用經驗不符。為解決該問題，后續需進一步收集系統使用初期和接近壽命極限這2 個時間階段的故障率數據。

隨著訓練樣本數量的增加，Kriging 模型的預測精度會逐漸提升，從已有的144 個樣本中隨機選取部分樣本作為訓練樣本訓練模型，并采用交叉驗證策略驗證模型預測精度。首先隨機選取20 個樣本作為訓練樣本訓練模型，并計算rRPRESS，隨后每次增加2 個訓練樣本數量訓練模型并計算rRPRESS，直到最終采用144個訓練樣本訓練模型，不同訓練樣本數量訓練模型對應的rRPRESS值如圖3 所示。從圖3 中可以看出，隨著訓練樣本數量增加，模型的rRPRESS值總體呈下降趨勢，即模型的預測精度逐漸提升，可說明后續若能進一步增加訓練樣本數量則能得到更加精確的預測結果。

圖3 訓練樣本數量對Kriging 故障率預測精度的影響Fig.3 Effect of training sample size on prediction accuracy of Kriging for failure rate

6 結語

本文對飛機系統的故障率敏感性分析問題進行了研究，結合Sobol 方法和Kriging 模型發展了一種適用于飛機系統故障率敏感性分析的方法。該方法通過收集不同使用年限和使用環境數據及對應的系統故障率，并將其作為訓練樣本訓練Kriging 模型成為滿足系統故障率預測要求的黑箱函數。然后通過利用訓練完成的Kriging 模型預測敏感性分析所需的樣本故障率，能有效避免真實故障率獲取困難的問題，因而可以更加高效地計算輸入變量的敏感性指標，并給出重要度排序，能很好地解決飛機系統故障率敏感性分析問題。最后，通過解決某飛機子系統的故障率敏感性分析問題，得到了影響系統故障率的各因素的主效應及全效應指標，證明了該方法在飛機系統故障率敏感性分析中的實際應用價值，可為飛機系統的外場維護使用及設計提升提供參考。