液體火箭發(fā)動(dòng)機(jī)自適應(yīng)克里金可靠性分析*

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

1. 浙江大學(xué) 流體動(dòng)力與機(jī)電系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，杭州 310027；2. 浙江大學(xué) 機(jī)械工程學(xué)院，杭州 310027

0 引言

近年來(lái)，隨著載人航天、探月工程等重大科技專項(xiàng)的實(shí)施，我國(guó)運(yùn)載火箭進(jìn)入了高強(qiáng)度發(fā)射期，對(duì)運(yùn)載火箭高可靠、高質(zhì)量等方面提出了更高的要求。液體火箭發(fā)動(dòng)機(jī)是運(yùn)載火箭等航天飛行器的心臟，需要具備更高的可靠度，因此對(duì)其運(yùn)行過(guò)程進(jìn)行故障監(jiān)測(cè)，從而在一定程度上提高運(yùn)載火箭的可靠性。液體火箭發(fā)動(dòng)機(jī)是一類高度復(fù)雜的動(dòng)力系統(tǒng)，由大量、相互之間存在廣泛聯(lián)系同時(shí)又進(jìn)行著不同工作過(guò)程的組件和部件構(gòu)成。隨著運(yùn)載火箭的有效載荷不斷增加，且外界工作環(huán)境存在不確定性，使發(fā)動(dòng)機(jī)工作在較為復(fù)雜的惡劣環(huán)境中。因此對(duì)其可靠性提出了更高的要求，液體火箭發(fā)動(dòng)機(jī)可靠性研究受到廣泛重視[1-3]。

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

基于上述結(jié)論，有學(xué)者提出了代理模型法，其主要原理為采用代理模型代替真實(shí)物理模型，運(yùn)用蒙特卡羅模擬計(jì)算失效概率。蒙特卡洛法需要較為龐大的樣本量，針對(duì)此，代理模型法提出將對(duì)應(yīng)樣本點(diǎn)輸入代理模型當(dāng)中，以較少的樣本數(shù)量進(jìn)行實(shí)時(shí)預(yù)測(cè)輸出，減少真實(shí)物理模型或耗時(shí)的有限元仿真模型的運(yùn)算次數(shù)與計(jì)算時(shí)間?，F(xiàn)有的代理模型主要包括：多項(xiàng)式響應(yīng)面模型[7]、神經(jīng)網(wǎng)絡(luò)模型[8]、支持向量機(jī)[9]以及克里金模型[10]。本文采用的代理模型為克里金模型，其優(yōu)勢(shì)在于其可以進(jìn)行全局代理，獲取了定義域內(nèi)樣本點(diǎn)區(qū)間的預(yù)測(cè)值，還包含預(yù)測(cè)局部方差，用于表征局部不確定性。

克里金模型對(duì)于模型的預(yù)測(cè)精確程度將直接影響分析結(jié)果的準(zhǔn)確性，因此基于較少的樣本構(gòu)建精確的克里金模型是當(dāng)前該領(lǐng)域的主要研究方向。Echard等[11]首先提出重要性采樣的思想，引入重要性密度分布函數(shù)提高采樣效率。在此基礎(chǔ)上，Xiao等[12]提出了分層重要性采樣的方法，進(jìn)一步提高采樣效率。此外Zhang等[13]提出的自適應(yīng)克里金面向重要度采樣方法和Wang等[14]提出的基于錯(cuò)誤率的自適應(yīng)克里金可靠性分析方法，為有效提高采樣效率，擬合精確的失效邊界提供新思路。但通過(guò)分析發(fā)現(xiàn)，現(xiàn)有方法大多存在樣本點(diǎn)聚集現(xiàn)象， 即同一區(qū)域內(nèi)樣本點(diǎn)數(shù)量過(guò)于密集，這將導(dǎo)致計(jì)算資源的浪費(fèi)。同時(shí)，現(xiàn)有方法的收斂條件不精準(zhǔn)，容易導(dǎo)致過(guò)采樣。

本文采用代理模型對(duì)液體火箭發(fā)動(dòng)機(jī)進(jìn)行可靠性分析。針對(duì)現(xiàn)有克里金代理模型存在的兩個(gè)弊端，提出了基于泰森多邊形區(qū)域劃分自適應(yīng)克里金的液體火箭發(fā)動(dòng)機(jī)可靠性分析方法。該方法能夠高效的構(gòu)建克里金代理模型，并用于預(yù)測(cè)液體火箭發(fā)動(dòng)機(jī)的失效概率。

1 基于代理模型的液體火箭發(fā)動(dòng)機(jī)可靠性分析

1.1 液體火箭發(fā)動(dòng)機(jī)可靠性分析現(xiàn)狀

液體火箭發(fā)動(dòng)機(jī)自身結(jié)構(gòu)極為復(fù)雜，不同組件之間的相互作用也會(huì)對(duì)運(yùn)行過(guò)程產(chǎn)生影響，存在隨機(jī)性。圖1給出了F-1發(fā)動(dòng)機(jī)的組成圖。針對(duì)液體火箭發(fā)動(dòng)機(jī)的可靠性分析，王婷婷等[1]提出威布爾分布和Bayes估計(jì)分析方法，研究液體火箭發(fā)動(dòng)機(jī)壽命分布參數(shù)以及該方法的可靠度。但這種方法未考慮液體火箭發(fā)動(dòng)機(jī)實(shí)際運(yùn)行時(shí)復(fù)雜環(huán)境給可靠性分析帶來(lái)的不確定性，無(wú)法得到實(shí)時(shí)的可靠性分析結(jié)果。

圖1 F-1液體火箭發(fā)動(dòng)機(jī)組成圖[15]

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

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

(1)

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

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

(2)

fcov(PF)為失效概率的差異系數(shù)，當(dāng)失效概率的差異系數(shù)值小于0.05時(shí)，可以認(rèn)為得到結(jié)果較為可信。

2 基于泰森多邊形區(qū)域劃分的自適應(yīng)克里金模型構(gòu)建方法

依據(jù)前文獲得的液體火箭發(fā)動(dòng)機(jī)在實(shí)際工作環(huán)境的失效概率，對(duì)其健康進(jìn)行監(jiān)測(cè)。計(jì)算失效概率計(jì)算失放概率的前提是構(gòu)建精確的液體火箭發(fā)動(dòng)機(jī)。針對(duì)傳統(tǒng)克里金代理模型構(gòu)建方法計(jì)算效率低，計(jì)算成本高等缺點(diǎn)，本文提出一種基于泰森多邊形區(qū)域劃分的自適應(yīng)克里金代理模型構(gòu)建方法，以較少的樣本構(gòu)建精確的克里金代理模型。

2.1 基于泰森多邊形的全局區(qū)域劃分

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

根據(jù)已有樣本點(diǎn)，生成所有相鄰兩個(gè)樣本點(diǎn)的垂直平分面，進(jìn)而由生成的垂直平分面將空間劃分成不同的區(qū)域，此區(qū)域?yàn)樘┥噙呅螀^(qū)域。圖2展示了一個(gè)由泰森多邊形劃分的二維案例。

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

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

2.2 改進(jìn)留一法

在全局區(qū)域劃分之后，需要選擇出在失效邊界附近的區(qū)域。該區(qū)域稱為最敏感區(qū)域，采用改進(jìn)留一法來(lái)實(shí)現(xiàn)。留一法[18](Leave-One-Out，簡(jiǎn)稱LOO)是一種模型驗(yàn)證優(yōu)化方法，常通過(guò)交叉驗(yàn)證檢驗(yàn)?zāi)Ｐ偷臏?zhǔn)確性。以下是留一法的具體流程：

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

(3)

將上述流程循環(huán)N次，最終得到所有樣本點(diǎn)的預(yù)測(cè)誤差。定義預(yù)測(cè)誤差最大的樣本點(diǎn)所屬的泰森多邊形區(qū)域?yàn)樽蠲舾袇^(qū)域。

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

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

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

2.3 主動(dòng)學(xué)習(xí)函數(shù)

在計(jì)算得到最敏感泰森多邊形區(qū)域之后，需在這個(gè)最敏感區(qū)域當(dāng)中選取下一個(gè)樣本點(diǎn)添加進(jìn)訓(xùn)練集，重新構(gòu)建克里金模型。本文采用的選擇樣本點(diǎn)方法為主動(dòng)學(xué)習(xí)函數(shù)。主動(dòng)學(xué)習(xí)函數(shù)用于確定構(gòu)建模型的下一個(gè)最佳采樣點(diǎn)，對(duì)于自適應(yīng)克里金模型應(yīng)用于可靠性分析領(lǐng)域必不可少。本文所采用的主動(dòng)學(xué)習(xí)函數(shù)是預(yù)期可行性函數(shù)(Expected Feasibility Function，EFF)。

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

(4)

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

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

2.4 整體流程

基于前文提及的方法，本文提出了一種基于泰森多邊形區(qū)域劃分自適應(yīng)克里金的液體火箭發(fā)動(dòng)機(jī)可靠性分析方法。該方法主要含有3個(gè)過(guò)程，包括全局泰森多邊形區(qū)域劃分，最敏感區(qū)域選擇和最有價(jià)值樣本點(diǎn)選擇。其具體流程如圖3所示。

圖3 基于泰森多邊形區(qū)域劃分自適應(yīng)克里金的液體火箭發(fā)動(dòng)機(jī)可靠性分析方法流程圖

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

步驟2：將初始樣本點(diǎn)X0及其對(duì)應(yīng)的輸出值作為訓(xùn)練集X，生成初始液體火箭發(fā)動(dòng)機(jī)克里金模型。具體包括：根據(jù)訓(xùn)練集樣本X，輸入到實(shí)際系統(tǒng)或有限元仿真模型中獲得其對(duì)應(yīng)的輸出Y；然后根據(jù)案例類型，設(shè)置克里金模型核函數(shù)類型，設(shè)置初始參數(shù)，設(shè)置參數(shù)優(yōu)化范圍；將訓(xùn)練集輸入，采用MATLAB的DACE工具箱構(gòu)建克里金模型；

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

步驟4：在全局內(nèi)進(jìn)行區(qū)域劃分，分為N個(gè)區(qū)域(N為現(xiàn)有訓(xùn)練集內(nèi)樣本數(shù)量)。全局區(qū)域劃分時(shí)采用泰森多邊形劃分策略，使得每個(gè)多邊形空間內(nèi)僅包含一個(gè)樣本點(diǎn)，且多邊形內(nèi)的任意位置離該多邊形的樣本點(diǎn)距離最近；

步驟5：根據(jù)改進(jìn)留一法，計(jì)算失效概率偏差，選出最敏感的區(qū)域C；

步驟6：從候選樣本集當(dāng)中劃分出屬于區(qū)域C內(nèi)的樣本點(diǎn)；

步驟7：采用主動(dòng)學(xué)習(xí)函數(shù)從區(qū)域C中選出一個(gè)樣本點(diǎn)，并加入訓(xùn)練集X，主動(dòng)學(xué)習(xí)函數(shù)采用的是預(yù)期可行性函數(shù)；

步驟8：根據(jù)現(xiàn)有訓(xùn)練集重新生成液體火箭發(fā)動(dòng)機(jī)克里金模型，具體細(xì)節(jié)如步驟2；

步驟9：判斷是否達(dá)到主動(dòng)學(xué)習(xí)函數(shù)收斂指標(biāo)，若未滿足則重復(fù)步驟4至步驟8。收斂指標(biāo)指的是式(5)，當(dāng)最敏感區(qū)域內(nèi)所有樣本點(diǎn)都滿足條件時(shí)認(rèn)為收斂。

max(fEF(x))≤0.001

(5)

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

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

步驟12：基于現(xiàn)有液體火箭發(fā)動(dòng)機(jī)克里金模型，根據(jù)蒙特卡羅模擬法計(jì)算系統(tǒng)失效概率Pf。

3 案例分析

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

(6)

式中，5個(gè)輸入變量分別是氧泵入口壓力X1、氧泵后溫度X2、氫泵入口壓力X3、氫泵后溫度X4以及燃?xì)獍l(fā)生器室壓X5均服從正態(tài)分布，均值0.5，方差0.1，5個(gè)變量相互獨(dú)立。

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

表1 測(cè)試案例結(jié)果對(duì)比

4 結(jié)論

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

后續(xù)工作可以從液體火箭發(fā)動(dòng)機(jī)時(shí)變可靠性分析研究著手，將外界載荷的時(shí)變特性納入考量，研究一種高效的時(shí)變可靠性分析方法。