液體火箭發動機自適應克里金可靠性分析*

胡偉飛 鄢繼銓 施欽楊 方健豪 劉振宇 譚建榮

1. 浙江大學 流體動力與機電系統國家重點實驗室，杭州 310027；2. 浙江大學 機械工程學院，杭州 310027

0 引言

近年來，隨著載人航天、探月工程等重大科技專項的實施，我國運載火箭進入了高強度發射期，對運載火箭高可靠、高質量等方面提出了更高的要求。液體火箭發動機是運載火箭等航天飛行器的心臟，需要具備更高的可靠度，因此對其運行過程進行故障監測，從而在一定程度上提高運載火箭的可靠性。液體火箭發動機是一類高度復雜的動力系統，由大量、相互之間存在廣泛聯系同時又進行著不同工作過程的組件和部件構成。隨著運載火箭的有效載荷不斷增加，且外界工作環境存在不確定性，使發動機工作在較為復雜的惡劣環境中。因此對其可靠性提出了更高的要求，液體火箭發動機可靠性研究受到廣泛重視[1-3]。

傳統可靠性分析的方法有：近似解析法、數值模擬法以及代理模型法，均適用于液體火箭發動機可靠性分析。近似解析法主要指的是一階可靠度法(FORM)[4]、二階可靠度法(SORM)[5]，其通過將設計點的失效域邊界線性化得到失效概率的估計。該類方法的缺點有：需要顯性液體火箭發動機數學模型，而且在非線性問題上應用效果較差。數值模擬法主要指蒙特卡羅模擬法[6]，該方法采用一定數量的樣本并獲取對應的響應值，基于這些樣本計算失效概率。計算結果的精度隨采樣的樣本點增加而提高。現有方法所需的計算成本和時間成本非常龐大。尤其針對液體火箭發動機這樣的關鍵裝備，獲取失效樣本的成本過高，收集大量數據不切實際。

基于上述結論，有學者提出了代理模型法，其主要原理為采用代理模型代替真實物理模型，運用蒙特卡羅模擬計算失效概率。蒙特卡洛法需要較為龐大的樣本量，針對此，代理模型法提出將對應樣本點輸入代理模型當中，以較少的樣本數量進行實時預測輸出，減少真實物理模型或耗時的有限元仿真模型的運算次數與計算時間。現有的代理模型主要包括：多項式響應面模型[7]、神經網絡模型[8]、支持向量機[9]以及克里金模型[10]。本文采用的代理模型為克里金模型，其優勢在于其可以進行全局代理，獲取了定義域內樣本點區間的預測值，還包含預測局部方差，用于表征局部不確定性。

克里金模型對于模型的預測精確程度將直接影響分析結果的準確性，因此基于較少的樣本構建精確的克里金模型是當前該領域的主要研究方向。Echard等[11]首先提出重要性采樣的思想，引入重要性密度分布函數提高采樣效率。在此基礎上，Xiao等[12]提出了分層重要性采樣的方法，進一步提高采樣效率。此外Zhang等[13]提出的自適應克里金面向重要度采樣方法和Wang等[14]提出的基于錯誤率的自適應克里金可靠性分析方法，為有效提高采樣效率，擬合精確的失效邊界提供新思路。但通過分析發現，現有方法大多存在樣本點聚集現象， 即同一區域內樣本點數量過于密集，這將導致計算資源的浪費。同時，現有方法的收斂條件不精準，容易導致過采樣。

本文采用代理模型對液體火箭發動機進行可靠性分析。針對現有克里金代理模型存在的兩個弊端，提出了基于泰森多邊形區域劃分自適應克里金的液體火箭發動機可靠性分析方法。該方法能夠高效的構建克里金代理模型，并用于預測液體火箭發動機的失效概率。

1 基于代理模型的液體火箭發動機可靠性分析

1.1 液體火箭發動機可靠性分析現狀

液體火箭發動機自身結構極為復雜，不同組件之間的相互作用也會對運行過程產生影響，存在隨機性。圖1給出了F-1發動機的組成圖。針對液體火箭發動機的可靠性分析，王婷婷等[1]提出威布爾分布和Bayes估計分析方法，研究液體火箭發動機壽命分布參數以及該方法的可靠度。但這種方法未考慮液體火箭發動機實際運行時復雜環境給可靠性分析帶來的不確定性，無法得到實時的可靠性分析結果。

圖1 F-1液體火箭發動機組成圖[15]

1.2 基于代理模型的可靠性分析

本文提出采用克里金模型法對液體火箭發動機進行可靠性分析。通過建立精確的液體火箭發動機克里金模型，預測液體火箭發動機的失效概率，采用的是傳統蒙特卡羅模擬方法。該方法通過采樣足夠數量的樣本，統計處于失效域的樣本數量，依照式(1)[16]計算出系統的失效概率PF。

(1)

式中，NG(x)<0定義為失效樣本數量；G代表構建的克里金代理模型；NMCS為采用蒙特卡洛法采樣樣本點總數量。

蒙特卡羅方法計算過程中所采用的樣本點越多，得到的結果越近似于真實值。在蒙特卡羅模擬計算中，得到的結果與樣本點數需要滿足如式(2)[16]所示的關系。

(2)

fcov(PF)為失效概率的差異系數，當失效概率的差異系數值小于0.05時，可以認為得到結果較為可信。

2 基于泰森多邊形區域劃分的自適應克里金模型構建方法

依據前文獲得的液體火箭發動機在實際工作環境的失效概率，對其健康進行監測。計算失效概率計算失放概率的前提是構建精確的液體火箭發動機。針對傳統克里金代理模型構建方法計算效率低，計算成本高等缺點，本文提出一種基于泰森多邊形區域劃分的自適應克里金代理模型構建方法，以較少的樣本構建精確的克里金代理模型。

2.1 基于泰森多邊形的全局區域劃分

本文所提出的方法采用了全局區域劃分策略，根據泰森多邊形的區域分割方式，進行空間區域劃分。具體流程如下：

根據已有樣本點，生成所有相鄰兩個樣本點的垂直平分面，進而由生成的垂直平分面將空間劃分成不同的區域，此區域為泰森多邊形區域。圖2展示了一個由泰森多邊形劃分的二維案例。

圖2 一個由泰森多邊形劃分的二維案例[17]

由于高維樣本點進行泰森多邊形分割后的多邊形單元是高維不規則邊界的多面體，因此很難對其進行數學描述。因此，本文采用蒙特卡羅方法生成大量隨機點從而近似描述高維問題中泰森多邊形區域。通過計算候選樣本點到空間中所有樣本點的距離，將每個候選樣本點劃分至已有的樣本點中最近的點。

2.2 改進留一法

在全局區域劃分之后，需要選擇出在失效邊界附近的區域。該區域稱為最敏感區域，采用改進留一法來實現。留一法[18](Leave-One-Out，簡稱LOO)是一種模型驗證優化方法，常通過交叉驗證檢驗模型的準確性。以下是留一法的具體流程：

假設空間中有N個樣本，每次剔除1個樣本，利用剩余的N-1個樣本重新構建模型，并采用被剔除的樣本作為測試集計算模型誤差，進行N次流程循環，具體如式(3)所示。

(3)

將上述流程循環N次，最終得到所有樣本點的預測誤差。定義預測誤差最大的樣本點所屬的泰森多邊形區域為最敏感區域。

以上是傳統留一法的過程。由于傳統留一法是將剩下的一個點作為測試集，通過比較在這個點的響應值的偏差來比較泰森多邊形區域的好壞，整個流程計算時間較長，計算成本較高。對于可靠性分析而言，只需要在失效邊界附近構建精確，能夠較好的擬合出極限狀態函數等于0(G(x)= 0)對應的曲面即可。因此不必構建一個全局精確的克里金模型，減少不必要的計算量和成本消耗。在失效邊界以外的區域，它的響應值不一定精確，只需要符號正確即可。

基于以上推論，本文所提出的方法對留一法進行修改，使其更適用于可靠性分析領域，具體流程如下：

假設空間中有N個樣本，每次剔除1個樣本，利用剩余的N-1個樣本構建克里金模型Gi。根據構建的克里金模型Gi計算剔除后的失效概率Pfi，共循環N次。除此之外采用N個樣本構建一個克里金模型，計算實際失效概率Pf。按照如式(3)，計算針對每個模型的失效概率偏差，最后將偏差最大模型對應的泰森多邊形區域選為最敏感區域。

2.3 主動學習函數

在計算得到最敏感泰森多邊形區域之后，需在這個最敏感區域當中選取下一個樣本點添加進訓練集，重新構建克里金模型。本文采用的選擇樣本點方法為主動學習函數。主動學習函數用于確定構建模型的下一個最佳采樣點，對于自適應克里金模型應用于可靠性分析領域必不可少。本文所采用的主動學習函數是預期可行性函數(Expected Feasibility Function，EFF)。

預期可行性函數fEF是由預期改進函數(Expected Improvement Function，EIF)改進而來的，由Jones等[19]在高效全局優化方法中首次提出。該預期改進函數是用來選擇一個新訓練樣本點的位置。但是預期改進函數更多應用于構建全局精確的模型。因此本文引入了預期可行性函數，以表示預期響應的真實值滿足等式約束G(x)=a的程度。其具體計算如式(4)[20]所示。

(4)

從敏感區域樣本點中選取預期可行性函數值最大的點作為下一個訓練樣本點。通過添加新的訓練樣本點，建立新的克里金代理模型，采用該模型計算另一組預期可行性函，進而由該預期可行性函數選擇另一個新的訓練點。以此類推，直到最大點的預期可行性函數值低于指定的公差。

失效概率估計的準確性總是隨著預期可行性函數停止準則的閾值的降低而增加的，但是代價是所需樣本點更多，增加了計算成本。本文采用最大的預期可行性函數值小于0.001作為迭代更新的收斂判據。

2.4 整體流程

基于前文提及的方法，本文提出了一種基于泰森多邊形區域劃分自適應克里金的液體火箭發動機可靠性分析方法。該方法主要含有3個過程，包括全局泰森多邊形區域劃分，最敏感區域選擇和最有價值樣本點選擇。其具體流程如圖3所示。

圖3 基于泰森多邊形區域劃分自適應克里金的液體火箭發動機可靠性分析方法流程圖

步驟1：根據不同輸入變量分布生成少量的初始樣本點X0，其中是采用拉丁超立方采樣(Latin Hypercube Sampling，LHS)方法生成初始樣本點；

步驟2：將初始樣本點X0及其對應的輸出值作為訓練集X，生成初始液體火箭發動機克里金模型。具體包括：根據訓練集樣本X，輸入到實際系統或有限元仿真模型中獲得其對應的輸出Y；然后根據案例類型，設置克里金模型核函數類型，設置初始參數，設置參數優化范圍；將訓練集輸入，采用MATLAB的DACE工具箱構建克里金模型；

步驟3：根據不同輸入變量分布生成一定數量候選樣本集S，生成候選樣本點使用的是拉丁超立方采樣方法；

步驟4：在全局內進行區域劃分，分為N個區域(N為現有訓練集內樣本數量)。全局區域劃分時采用泰森多邊形劃分策略，使得每個多邊形空間內僅包含一個樣本點，且多邊形內的任意位置離該多邊形的樣本點距離最近；

步驟5：根據改進留一法，計算失效概率偏差，選出最敏感的區域C；

步驟6：從候選樣本集當中劃分出屬于區域C內的樣本點；

步驟7：采用主動學習函數從區域C中選出一個樣本點，并加入訓練集X，主動學習函數采用的是預期可行性函數；

步驟8：根據現有訓練集重新生成液體火箭發動機克里金模型，具體細節如步驟2；

步驟9：判斷是否達到主動學習函數收斂指標，若未滿足則重復步驟4至步驟8。收斂指標指的是式(5)，當最敏感區域內所有樣本點都滿足條件時認為收斂。

max(fEF(x))≤0.001

(5)

步驟10：計算當前克里金模型下的失效概率，計算方法采用的是蒙特卡羅模擬法；

步驟11：根據式(2)判斷是否滿足蒙特卡羅法收斂指標失效概率協方差是否小于0.05，若未滿足，則根據生成新的候選樣本加入候選樣本集S，并返回步驟3；

步驟12：基于現有液體火箭發動機克里金模型，根據蒙特卡羅模擬法計算系統失效概率Pf。

3 案例分析

本章采用的測試案例是基于楊晉朝[21]提供的液體火箭發動機失效數據，得到的液體火箭發動機性能方程，如式(6)所示：

(6)

式中，5個輸入變量分別是氧泵入口壓力X1、氧泵后溫度X2、氫泵入口壓力X3、氫泵后溫度X4以及燃氣發生器室壓X5均服從正態分布，均值0.5，方差0.1，5個變量相互獨立。

在測試過程中，初始樣本點選擇為12個，每次生成的候選樣本集數量為105個。將本文所提出的方法與傳統的自適應克里金模型-蒙特卡羅模擬方法[20](An active learning reliability method combining Kriging and Monte Carlo Simulation, AK-MCS)作比較，并采用蒙特卡羅模擬結果作為真實值。并通過5次仿真測試后求取平均結果，驗證本方法的魯棒性。最終的實驗對比結果如表1所示。表中展示了兩種方法所需的樣本點個數和誤差率。從表中可以看出，本文所提出方法與傳統自適應采樣方法的計算誤差都小于1%，達到精度計算要求。此外，本文提出方法所需的樣本點個數僅為163.4個，而AK-MCS方法需要314.2個樣本點。在本案例中，本文所提出方法效率比AK-MCS方法提高了48%，相對于傳統方法能夠快速、準確地對液體火箭發動機進行可靠性分析。

表1 測試案例結果對比

4 結論

針對傳統可靠性分析方法計算效率低的問題，本文提出了一種基于泰森多邊形區域劃分自適應克里金的液體火箭發動機可靠性分析方法。利用泰森多邊形進行全局區域劃分，結合改進留一法對區域進行篩選，極大程度的提高了模型構建的效率，能夠用較少的樣本點構建同樣精度要求的液體火箭發動機模型。同時，該模型能夠結合外界載荷的不確定性，對發動機進行健康監測和故障檢測。在理論測試中，基于本文提出方法構建的代理模型測試計算誤差與其他模型相比小于1%，效率提高了48%。本文側重于理論方法研究，如有真實液體火箭發動機模型或有限元仿真模型，采用本方法也能對其進行可靠性分析，計算出特定工作狀況下的失效概率。

后續工作可以從液體火箭發動機時變可靠性分析研究著手，將外界載荷的時變特性納入考量，研究一種高效的時變可靠性分析方法。