張海瑞，王 浩，王 堯，洪東跑

(中國運載火箭技術(shù)研究院，北京 100076)

0 引 言

飛行器分離是飛行器發(fā)射任務(wù)的關(guān)鍵環(huán)節(jié)之一，其功能是在一定的高度將完成預(yù)定功能的部分結(jié)構(gòu)分離并拋去，旨在改善總體性能，增大飛行速度，提高運載能力[1-3]。飛行器實施分離時的惡劣環(huán)境將導(dǎo)致兩體可靠分離困難，嚴(yán)重制約飛行器可靠性水平的提升。在分離過程中存在著諸多不確定性因素[4]，如反推發(fā)動機作用時間不同步，下面級殘余推力偏斜以及殘余推力作用點位置偏差等。在這些不確定性因素的綜合影響下，兩體分離存在碰撞風(fēng)險，可能導(dǎo)致分離任務(wù)失敗。

傳統(tǒng)的分離分析方法通常是基于參數(shù)偏差進行的，通過偏差組合考慮最嚴(yán)酷的情況[5]，以此判斷分離設(shè)計方案是否滿足要求。然而隨著飛行器可靠性要求不斷提高以及新技術(shù)方法的發(fā)展應(yīng)用，傳統(tǒng)方法為了包絡(luò)最嚴(yán)酷的極限條件，其設(shè)計方案往往偏于保守，增加結(jié)構(gòu)重量，在某些情況下甚至嚴(yán)重影響飛行器總體性能，難以滿足精細(xì)化設(shè)計需求。

針對這一問題，目前的研究趨勢是綜合考慮不確定性因素影響，將不確定性注入分離動力學(xué)仿真模型，建立飛行器分離可靠性模型，實現(xiàn)分離方案精細(xì)化分析[6-8]。然而，不確定性因素和模型非線性的影響導(dǎo)致分離可靠性求解困難。國內(nèi)外學(xué)者系統(tǒng)研究了利用蒙特卡洛方法分析不確定性因素對分離運動的影響[9-10]，Roshanian等[11]采用蒙特卡洛方法研究了不確定性因素對分離運動的影響，給出了分離運動參數(shù)的包絡(luò)曲線。沙建科等[12]針對某型導(dǎo)彈級間冷分離方案，采用蒙特卡洛方法對隨機不確定性進行了分離仿真分析，給出了兩體不發(fā)生碰撞的分離可靠性。蒙特卡洛方法具有無偏性以及非侵入性，在工程中有著廣泛的應(yīng)用。當(dāng)前，對如何進行分離可靠性高效分析與精確評估的研究較少，特別是隨著飛行器分離可靠性指標(biāo)的不斷提高，蒙特卡洛方法調(diào)用分離動力學(xué)仿真模型的次數(shù)也會隨之急劇增加，從而消耗巨大的計算資源，甚至在有限的研制時間內(nèi)無法給出精確的分離可靠性分析結(jié)果。

為進一步提高可靠性分析的效率，國內(nèi)外學(xué)者提出了諸多高效的可靠性分析方法，其中基于代理模型的可靠性分析方法具有簡單直接，適用廣泛的特點，是可靠性分析研究的發(fā)展趨勢之一[7]。文獻[13]基于Kriging模型的特點提出了高效全局優(yōu)化(Efficient global optimization, EGO)方法用于求解復(fù)雜函數(shù)的優(yōu)化問題。文獻[14]基于EGO方法提出了高效全局可靠性分析(Efficient global reliability analysis, EGRA)方法求解復(fù)雜問題的可靠度，其中，采用期望可行函數(shù)(Expected feasibility function,EFF)序貫增加訓(xùn)練點。文獻[15]結(jié)合Kriging模型的特點通過U函數(shù)序貫增加訓(xùn)練點，提出了主動學(xué)習(xí)Kriging方法求解隱式功能函數(shù)的可靠度。文獻[16]提出了風(fēng)險期望函數(shù)(Expected risk function，ERF)序貫加點，求解非線性功能函數(shù)的混合可靠度。文獻[17]針對小失效概率問題，將主動學(xué)習(xí)Kriging和子集模擬方法相結(jié)合，提出了AK-SS方法求解隱式功能函數(shù)的小失效概率問題。目前，可靠性分析方法已經(jīng)取得了長足的進步，但是面向高超聲速飛行器分離任務(wù)需求，開展高效分離可靠性分析的研究還較少。

本文以某軸對稱式飛行器級間冷分離方案為研究對象，建立了分離動力學(xué)仿真模型，為滿足精細(xì)化設(shè)計需求，通過綜合考慮分離過程不確定性因素的影響，構(gòu)建了飛行器分離可靠性模型。結(jié)合代理模型和蒙特卡洛方法的特點，提出了改進主動學(xué)習(xí)Kriging(Improved active kriging, IAK)的分離可靠性分析方法，給出了一種基于不確定性的飛行器分離可靠性建模與分析方法。

1 飛行器分離動力學(xué)建模

發(fā)射過程中飛行器的正常分離通常包括助推器分離、級間分離、整流罩分離、星箭分離等[1]，在飛行器設(shè)計過程中，級間分離技術(shù)仍然是需要攻克的一項關(guān)鍵技術(shù)，包括級間熱分離和級間冷分離兩種方法[18]。本文選取典型的軸對稱式飛行器級間冷分離方案為研究對象。針對其環(huán)境特點與性能要求，建立分離動力學(xué)模型。分離過程中，上面級受到重力、分插拔脫力和氣動力的影響，下面級受到重力、分插拔脫力、主發(fā)動機殘余推力、反推發(fā)動機推力以及氣動力的影響，飛行器級間冷分離示意圖如圖1所示。由于上面級發(fā)動機位于下面級的殼體內(nèi)，故在分離過程中需要著重關(guān)注兩者之間的相對距離。

圖1 飛行器級間分離示意圖Fig.1 Illustration of vehicle stage separation

定義飛行器上面級和下面級的彈體坐標(biāo)系原點與各自質(zhì)心重合，x軸與各自彈體縱軸重合；定義分離坐標(biāo)系與分離初始時刻飛行器組合體的彈體坐標(biāo)系重合，由于忽略地球自轉(zhuǎn)影響，分離坐標(biāo)系即為分離慣性參考系。以此為基礎(chǔ)，在分離坐標(biāo)系下分別針對上面級和下面級建立質(zhì)心運動方程如下

(1)

其中，m是分析對象的質(zhì)量；vx,vy,vz分別是分析對象的速度矢量在分離坐標(biāo)系下的速度投影；Fx,Fy,Fz分別是分析對象受到的合力在分離坐標(biāo)系下的投影。進一步給出分析對象彈體坐標(biāo)系到分離坐標(biāo)系的轉(zhuǎn)換矩陣為

GB=

(2)

其中，φ,ψ,γ分別為分析對象的彈體坐標(biāo)系相對分離坐標(biāo)系的歐拉角，通過上述轉(zhuǎn)換矩陣將分析對象在彈體坐標(biāo)系下受到的合力轉(zhuǎn)換為分離坐標(biāo)系下受到的合力。

進一步，分別在上面級和下面級的彈體坐標(biāo)系下建立繞質(zhì)心的運動方程為

(3)

其中，ωx1,ωy1,ωz1為分析對象轉(zhuǎn)動角速度在彈體坐標(biāo)系中的分量；Ix1,Iy1,Iz1為分析對象相對其彈體坐標(biāo)系的轉(zhuǎn)動慣量。分別以上面級和下面級為分析對象，補充建立速度與位置、歐拉角與角速度之間的關(guān)系，如式(4)所示:

(4)

針對分離過程的具體特點，結(jié)合上述方程分別建立組合體、上面級和下面級的六自由度剛體運動及動力學(xué)仿真模型?？紤]到下面級主發(fā)動機關(guān)機后仍會產(chǎn)生殘余推力，且殘余推力隨時間變化而逐漸減小，進一步考慮殘余推力的偏移偏斜影響，如圖2所示。

P是殘余推力的等效作用點，X10,Y10,Z10與下面級彈體坐標(biāo)系方向平行。其中，λ角為殘余推力的推力線與下面級彈體中軸線的夾角；θ角為殘余推力的推力線在截面Y10PZ10上的投影與Y10軸的夾角，給出殘余推力在彈體坐標(biāo)系下的分量為

(5)

圖2 殘余推力偏移偏斜示意圖Fig.2 Illustration of the offset and deviation of residual thrust

由于偏移偏斜的影響，殘余推力的作用點與下面級質(zhì)心的相對位置為rpc=[xpc,ypc,zpc]T，發(fā)動機推力偏移偏斜產(chǎn)生的力矩為

(6)

在飛行器分離過程中，存在多種分離故障失效模式，其中較為典型的故障判據(jù)為兩體在分離時發(fā)生干涉碰撞，這里采用文獻[10]中的碰撞檢測方案計算上面級噴管底端面中心相對下面級中軸線的最大偏移距離d，根據(jù)分離結(jié)構(gòu)的具體特點，給定最大偏移距離的閾值d0，若從開始分離到噴管完全脫離底端面的過程中，最大偏移距離小于給定閾值d0，則判定分離過程中未發(fā)生碰撞。中心點偏移示意圖如圖3所示。

圖3 中心點偏移示意圖Fig.3 Illustration of the center offset

(7)

2 飛行器分離可靠性建模

結(jié)合分離動力學(xué)模型，針對軸對稱式飛行器級間冷分離方案的特點，綜合考慮質(zhì)量、氣動、殘余推力以及反推發(fā)動機等不確定性因素影響，給出分離過程的不確定性因素如表1所示。

表1 不確定性因素名稱及符號Table 1 Names and symbols of uncertainty parameters

通過中心差分方法對上述20個分離不確定性參數(shù)進行顯著性分析[19-20]，經(jīng)歸一化處理后，給出上面級噴管底端面中心相對下面級中軸線最大偏移距離d影響最大的10個因素如圖4所示。

圖4 不確定性因素對最大偏移距離影響百分比排列圖Fig.4 Effects of uncertainty parameters on maximum offset distance in the separation process

基于上述分析，選取主要不確定性因素包括：反推發(fā)動機不同步時間t0(x16)，氣動特性偏差Wqd(x20)，反推推力作用點位置xFt(x14)，殘余推力的推力線與下面級彈體中軸線的夾角λ(x18)以及殘余推力作用點位置xFe(x17)。這5個主要不確定性因素的影響百分比總和超過95%，識別主要不確定性因素，濾除影響較小的次要不確定性因素。

為滿足精細(xì)化設(shè)計需求，在建立分離不確定性模型的基礎(chǔ)上，將主要不確定性因素注入分離動力學(xué)模型，建立飛行器分離可靠性模型。結(jié)合分離碰撞判據(jù)，給出可靠性模型的極限狀態(tài)方程為分離過程上面級噴管底端面的最大偏移距離d與給定閾值d0之差等于零，即：

g(t0,xFt,xFe,λ,Wqd)=d(t0,xFt,xFe,λ,Wqd)-

d0=0

(8)

其中，t0,xFt,xFe,λ,Wqd為不確定性參數(shù)，失效域定義為

D={(t0,xFt,xFe,λ,Wqd)g(t0,xFt,xFe,λ,Wqd)≥0}

(9)

因而在分離過程中兩體不發(fā)生碰撞的分離可靠性為

R=Pr(d(t0,xFt,xFe,λ,Wqd)-d0<0)

(10)

進一步定義可靠性模型極限狀態(tài)方程的聯(lián)合概率密度函數(shù)為p(t0,xFt,xFe,λ,Wqd)，分離可靠性還可以表示為

(11)

3 飛行器分離可靠性分析

由于不確定性參數(shù)為隨機變量，通過試驗數(shù)據(jù)和工程經(jīng)驗確定分離不確定參數(shù)的分布，進而由式(11)可得分離可靠度。鑒于不確定參數(shù)較多，難以獲得失效域D與聯(lián)合概率密度函數(shù)p(t0,xFt,xFe,λ,Wqd)的顯示表達(dá)，直接積分求解難度較大且復(fù)雜耗時，因而通常利用蒙特卡洛法進行可靠性分析，此方法精度較高，但隨著分離可靠性指標(biāo)的提高，調(diào)用分離仿真模型的次數(shù)也會隨之急劇增加，導(dǎo)致求解效率較低，難以滿足快速設(shè)計迭代的工程需求。針對這一問題，文獻[15,21]提出了主動學(xué)習(xí)Kriging的可靠性分析方法，具有很高的效率，在蒙特卡洛樣本點中直接選取失效概率最大的樣本點作為新增訓(xùn)練點，然而其新增訓(xùn)練點的選擇完全取決于蒙特卡洛樣本點集中失效概率最大的樣本點，無法選取樣本點集之外失效概率更大的未知點，從而限制了效率和精度的進一步提升。

結(jié)合代理模型和蒙特卡洛方法的特點，提出了一種改進主動學(xué)習(xí)Kriging的分離可靠性分析方法，通過新的采樣策略在失效概率最大樣本點附近完成采樣，在滿足采樣點與訓(xùn)練點集相關(guān)條件的基礎(chǔ)上，給出失效概率更大的新增訓(xùn)練點，使得新增訓(xùn)練點的最大偏移距離更加趨于給定閾值。

記分離不確定性參數(shù)向量x=(t0,xFt,xFe,λ,Wqd)，利用Kriging代理模型，將真實最大偏移距離d與給定閾值d0之差表達(dá)為如下形式

(12)

其中，f(x)=[f1(x),f2(x),…,fn(x)]T為回歸基函數(shù)，針對分離動力學(xué)模型的具體特點，采用二階多項式回歸模型；β=[β1,β2,…,βn]T為回歸系數(shù)，z(x)為高斯隨機修正過程。在給出未知點x預(yù)測最大偏移距離與給定閾值之差的基礎(chǔ)上，還可以進一步給出預(yù)測值的標(biāo)準(zhǔn)差，其預(yù)測值服從如下高斯分布

(13)

其中，μG(x)為Kriging的預(yù)測值，σG(x)為預(yù)測值的標(biāo)準(zhǔn)差。

(14)

式中:Rt為目標(biāo)分離可靠性。在樣本點集XMC中采用如下學(xué)習(xí)函數(shù)確定失效概率最大的樣本點A。

(15)

其中，μ(x)為預(yù)測的分離最大偏移距離與給定閾值之差，σ(x)為預(yù)測值的標(biāo)準(zhǔn)差。當(dāng)學(xué)習(xí)函數(shù)的最小值大于2時，所有樣本點的最大失效概率僅為Φ(-2)≈0.0228，滿足分離可靠性的計算需求，即滿足收斂條件。

式(15)的U函數(shù)更傾向于選取最大偏移距離在閾值附近或者尚未探索區(qū)域的樣本點。因而，進一步將A點分為兩類：若A點預(yù)測最大偏移距離在給定閾值附近，記為A1；若A點預(yù)測最大偏移距離遠(yuǎn)離給定閾值，則屬于探索未知區(qū)域的情況，記為A2，如圖5所示。

圖5 采樣策略示意圖Fig.5 Illustration of the proposed sampling strategy

針對預(yù)測最大偏移距離在給定閾值附近的情況，采用新的采樣策略確定失效概率更大的新增訓(xùn)練點，采樣策略如下：首先通過A1點的符號判斷其屬于失效點還是可行點，若A1點屬于可行點，則選取失效點集為樣本點集XB，反之，選取可行點集為樣本點集XB。之后，在樣本點集XB中給出與A1點距離最近的樣本點B，如圖5所示。最后，連接A1,B兩點，取N+1個采樣點將線段AB等分為N段，記采樣點集為XP。在采樣點集XP中選取失效概率最大的采樣點作為新增訓(xùn)練點。

針對預(yù)測最大偏移距離遠(yuǎn)離給定閾值的情況，其目的是探索未知區(qū)域，若在給定閾值附近選取失效概率更大的采樣點，此時依然無法提高代理模型未知區(qū)域的精度，因而直接選取失效概率密度最大的樣本點A2為新增訓(xùn)練點。

由于采樣點集XP與當(dāng)前訓(xùn)練點集XT不相關(guān)，因而，存在迭代后期新增訓(xùn)練點xnew與訓(xùn)練點集XT相關(guān)程度過大而影響Kriging模型精度的情況，需要進一步補充相關(guān)條件。記當(dāng)前訓(xùn)練點集XT={x1,x2,…,xm}，則新增訓(xùn)練點xnew的Kriging修正函數(shù)表述如下

z(xnew)=rT(xnew)Q-1(y-F·β)

(16)

其中，y=[g(x1),g(x2),…,g(xm)]T為訓(xùn)練點集的最大偏移距離與給定閾值之差；Q為相關(guān)矩陣，rT(x)為新增訓(xùn)練點xnew與訓(xùn)練點集XT之間的相關(guān)向量，定義為

r(x)=[R(xnew,x1),R(xnew,x2),…,R(xnew,xm)]

(17)

其中,R為相關(guān)函數(shù)。補充兩者之間的相關(guān)條件如下

max{r(xnew,XT)}<ε

(18)

其中,ε為相關(guān)程度的閾值，通常ε∈[0.9, 1]。通過上述方法，將新增訓(xùn)練點與訓(xùn)練點集之間的擁擠距離轉(zhuǎn)化為兩者之間的相關(guān)函數(shù)，在不影響Kriging模型精度的基礎(chǔ)上，選取了失效概率更大的新增訓(xùn)練點，建立了分離可靠性極限狀態(tài)方程的局部高精度代理模型，結(jié)合蒙特卡洛方法，減少分離模型的調(diào)用次數(shù)，提高分離可靠性分析效率。

綜上所述，改進主動學(xué)習(xí)Kriging方法的流程圖如圖6所示。

圖6 改進主動學(xué)習(xí)Kriging流程圖Fig.6 Flowchart of the improved active learning Kriging approach

4 實例結(jié)果

以某低空高速軸對稱式飛行器級間冷分離方案為應(yīng)用對象，飛行高度40 km，分離沖量由兩個對稱安裝的反推火箭提供?？紤]加工及安裝誤差的影響，主要不確定性因素分布類型及分布參數(shù)如表2所示。

根據(jù)識別的分離不確定性模型隨機選取50000個蒙特卡洛樣本點，采用相同的樣本點集分別采用主動學(xué)習(xí)Kriging和改進主動學(xué)習(xí)Kriging進行級間分離可靠性分析。為進一步驗證改進主動學(xué)習(xí)Kriging的穩(wěn)健性，給出10次分離可靠性分析結(jié)果，如表3所示。

改進主動學(xué)習(xí)Kriging的分離仿真模型調(diào)用次數(shù)更少，分析效率更高，其可靠性分析效率平均提高了16.34%。隨機選取10次不同的初始訓(xùn)練點集和蒙特卡洛樣本點集，主動學(xué)習(xí)Kriging和改進主動學(xué)習(xí)Kriging結(jié)果的標(biāo)準(zhǔn)差近似相同，校驗了改進主動學(xué)習(xí)Kriging方法的穩(wěn)健性。給出其中第7次仿真訓(xùn)練點集的分離最大偏移距離與給定閾值之差的歷程對比以及U函數(shù)迭代歷程對比，如圖7、圖8所示。

表2 主要不確定性因素分布類型及其參數(shù)Table 2 Distribution models and parameters of the significant uncertainty parameters

主動學(xué)習(xí)Kriging(AK)和改進主動學(xué)習(xí)Kriging (IAK)都趨于選擇分離最大偏移距離在給定閾值附近的樣本點。因而，兩者在極限狀態(tài)方程附近都具有較高的精度。同時，U函數(shù)迭代歷程表明，改進主動學(xué)習(xí)Kriging方法提高了采樣策略的加點質(zhì)量，尤其是早期加點質(zhì)量，使得U函數(shù)更早地進入到上升階段，提高了可靠性分析效率。進一步，采用蒙特卡洛方法計算分離可靠性，驗證方法的準(zhǔn)確性和有效性，給出分離可靠性分析對比如表4所示。其中仿真平臺的處理器為Intel(R) Core(TM) i7-6700U CPU @3.40 GHz,內(nèi)存4.00 GB。

表3 飛行器級間分離可靠性分析結(jié)果Table 3 Vehicle stage separation reliability analysis results

表4 飛行器級間分離可靠性分析對比Table 4 Comparison of reliability analysis of vehicle stage separation

圖7 新增訓(xùn)練點最大偏移距離與給定閾值之差歷程對比Fig.7 Iteration history of the difference between maximum offset distance and given threshold for new training points

圖8 U函數(shù)迭代歷程對比Fig.8 Comparison of U function iterations

為進一步提升分離任務(wù)的可靠性，考慮不同高度對分離可靠性的影響，采用上述飛行器可靠性建模與分析方法給出33 km～44 km之間飛行器兩體不發(fā)生碰撞的分離可靠性，如圖9所示。

針對低空高速軸對稱式飛行器級間冷分離方案，飛行高度對分離可靠性具有顯著影響，由于空氣密度較大，飛行馬赫數(shù)較高，飛行動壓較大，因而氣動力及氣動力矩的影響顯著。選取適當(dāng)?shù)姆蛛x高度可以提高分離任務(wù)的可靠性，保證兩體可靠分離，在量化分離不確定性的同時，給出分離可靠性定量分析結(jié)果，為提升分離方案的可靠性提供了重要的方法支撐。

圖9 不同高度下的分離可靠性Fig.9 Separation reliability at different altitudes

5 結(jié) 論

1)本文以軸對稱式飛行器級間冷分離方案為研究對象，通過考慮分離過程中不確定性因素的影響，將主要不確定性注入分離動力學(xué)仿真模型,建立了飛行器分離可靠性模型，實現(xiàn)了分離任務(wù)可靠性的定量分析，提升了分離可靠性的精度。

2)針對飛行器分離可靠性模型的特點，提出了改進主動學(xué)習(xí)Kriging的分離可靠性分析方法，提升了分離可靠性定量分析的效率。

3)結(jié)合實例結(jié)果驗證了改進主動學(xué)習(xí)Kriging方法的正確性、高效性和穩(wěn)健性，并有效減少分離可靠性分析的計算時間。下一步的研究重點是與多點并行加點策略相結(jié)合，進一步降低計算時間，提高擬合效率。