祝學(xué)軍，卜奎晨，王 浩，高 峰，趙長(zhǎng)見(jiàn)

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

飛行安全區(qū)，又稱飛行安全控制區(qū)，是指火箭實(shí)際飛行時(shí)，可能使發(fā)射場(chǎng)、航跡區(qū)的人員、設(shè)施遭受損傷和破壞的區(qū)域[1]。根據(jù)火箭正常飛行或故障飛行情況，可將安全區(qū)劃分為發(fā)射場(chǎng)安全區(qū)(或稱為首區(qū)安全區(qū))、子級(jí)殘骸安全區(qū)、航區(qū)安全區(qū)及落區(qū)安全區(qū)。由于飛行安全區(qū)是火箭正常飛行時(shí)分離殘骸或者異常飛行時(shí)的自毀殘骸可能到達(dá)的區(qū)域，該散布區(qū)域大小與火箭性能和飛行環(huán)境密切相關(guān)，因此在進(jìn)行安全區(qū)計(jì)算時(shí)需要考慮各種偏差和干擾的影響，如發(fā)動(dòng)機(jī)性能偏差、結(jié)構(gòu)質(zhì)量偏差、氣動(dòng)偏差、大氣偏差及風(fēng)干擾等。

傳統(tǒng)彈道式飛行器采用拋物線彈道，大部分彈道位于真空中，飛行器在大氣層內(nèi)飛行的時(shí)間較短，氣動(dòng)偏差、大氣偏差及風(fēng)干擾等因素對(duì)安全區(qū)范圍影響相對(duì)較小，傳統(tǒng)彈道式飛行器在進(jìn)行飛行安全區(qū)計(jì)算時(shí)，往往采用幾種主要偏差極限疊加的方法進(jìn)行安全區(qū)計(jì)算，計(jì)算方法偏保守。由于我國(guó)經(jīng)濟(jì)不斷發(fā)展，人員經(jīng)營(yíng)和活動(dòng)范圍不斷擴(kuò)大，航天發(fā)射可用區(qū)域不斷縮小，導(dǎo)致發(fā)射場(chǎng)安全控制實(shí)施難度不斷加大。傳統(tǒng)極限偏差疊加的安全區(qū)計(jì)算方法因過(guò)于保守已不再適用。

相比傳統(tǒng)彈道式飛行器，助推滑翔飛行器的助推火箭殘骸在大氣層中飛行的時(shí)間更長(zhǎng)，受大氣環(huán)境和氣動(dòng)力偏差的影響更嚴(yán)重，其分離殘骸安全區(qū)范圍顯著增大，進(jìn)一步加大了發(fā)射場(chǎng)安全控制實(shí)施難度，此外，還存在對(duì)發(fā)射時(shí)間選擇等方面的需求，對(duì)安全區(qū)精細(xì)化設(shè)計(jì)需求更高。

國(guó)內(nèi)外學(xué)者在飛行安全區(qū)精細(xì)化計(jì)算方面開(kāi)展了較多的研究工作[2-4]，并取得了一定的成果。文獻(xiàn)[1]采用Monte Carlo方法進(jìn)行了運(yùn)載火箭殘骸落區(qū)的計(jì)算，有效地減少了安全防范區(qū)域。文獻(xiàn)[2]針對(duì)火箭在不同季節(jié)、不同月份發(fā)射的飛行包絡(luò)開(kāi)展了打靶仿真，用于確定發(fā)射時(shí)間。但由于偏差組合工況過(guò)多、仿真時(shí)間不充裕，無(wú)法滿足使用的需求，為此采取了選取惡劣工況的方法來(lái)降低仿真次數(shù)。可見(jiàn)，Monte Carlo方法將安全區(qū)范圍、偏差項(xiàng)及發(fā)生概率建立了關(guān)系，是一種較為精細(xì)化的安全區(qū)計(jì)算方法，但當(dāng)計(jì)算工況較多時(shí)，存在仿真時(shí)間過(guò)長(zhǎng)的問(wèn)題。國(guó)內(nèi)外學(xué)者針對(duì)如何確定Monte Carlo方法打靶次數(shù)開(kāi)展研究[5-6]，并得出“對(duì)于一般的工程技術(shù)問(wèn)題，打靶次數(shù)3000～5000即可滿足工程精度要求”的結(jié)論。但由于存在需要單獨(dú)施加偏差(如風(fēng)干擾)及需對(duì)發(fā)射時(shí)間(不同季節(jié)、不同月份)進(jìn)行選擇等問(wèn)題，打靶次數(shù)仍然會(huì)大量增加。為了進(jìn)一步加快仿真速度，近年來(lái)在采用代理模型加速M(fèi)onte Carlo打靶方面也有一定的研究，文獻(xiàn)[7-8]分別針對(duì)Kriging模型開(kāi)展了相關(guān)研究，并取得了一定的成果。

本文結(jié)合Monte Carlo方法和代理模型的特點(diǎn)，提出了一種基于優(yōu)化加點(diǎn)Kriging模型的助推火箭殘骸安全區(qū)分析方法，為助推火箭殘骸等安全區(qū)預(yù)示提供理論參考。

1 彈道計(jì)算模型

1.1 運(yùn)動(dòng)方程

火箭及分離殘骸的彈道計(jì)算質(zhì)心運(yùn)動(dòng)動(dòng)力學(xué)方程和繞質(zhì)心轉(zhuǎn)動(dòng)動(dòng)力學(xué)方程[9]見(jiàn)式(1)和式(2)。

(1)

(2)

質(zhì)心運(yùn)動(dòng)學(xué)方程見(jiàn)式(3)。

(3)

式中，x、y、z為發(fā)射系位置分量，Vx、Vy、Vz為發(fā)射系速度分量。

為防止姿態(tài)角解算出現(xiàn)奇異問(wèn)題，采用四元數(shù)方法進(jìn)行姿態(tài)角解算處理。基于四元數(shù)的繞質(zhì)心轉(zhuǎn)動(dòng)運(yùn)動(dòng)學(xué)方程見(jiàn)式(4)和式(5)。

(4)

(5)

式中，q0、q1、q2、q3為四元數(shù)，φT、ψT、γT為姿態(tài)角。

1.2 主要偏差情況

在針對(duì)助推火箭殘骸散布區(qū)域等進(jìn)行安全區(qū)計(jì)算時(shí)，需要考慮助推段發(fā)動(dòng)機(jī)性能偏差、質(zhì)量特性偏差、氣動(dòng)系數(shù)偏差、大氣偏差(溫度、密度、壓強(qiáng))、風(fēng)干擾等不確定干擾因素，同時(shí)考慮殘骸飛行段的氣動(dòng)系數(shù)偏差、大氣偏差及風(fēng)干擾等不確定干擾因素。干擾因素的分布規(guī)律依賴于統(tǒng)計(jì)數(shù)據(jù)，一般可以認(rèn)為發(fā)動(dòng)機(jī)性能偏差、質(zhì)量特性偏差、氣動(dòng)系數(shù)偏差、大氣偏差等服從高斯分布，記為：

(6)

發(fā)動(dòng)機(jī)性能偏差一般包括發(fā)動(dòng)機(jī)推力偏斜、發(fā)動(dòng)機(jī)工作時(shí)間偏差、推力偏差和總沖偏差等。

風(fēng)干擾一般取為常值項(xiàng)，風(fēng)速大小VW隨飛行高度變化，計(jì)算公式見(jiàn)式(7)。由于風(fēng)向或射向的不確定性，通常按照順風(fēng)、逆風(fēng)、左側(cè)風(fēng)、右側(cè)風(fēng)四個(gè)工況分別施加并針對(duì)其他偏差項(xiàng)進(jìn)行模擬打靶分析[3]。

VW=f(h)

(7)

2 安全區(qū)預(yù)示方法

2.1 安全區(qū)預(yù)示模型

在不確定性干擾因素的影響下，工程上殘骸安全區(qū)范圍通常假設(shè)為某一矩形區(qū)域，由縱向最遠(yuǎn)距離εymax、縱向最近距離εymin、橫向最遠(yuǎn)距離εxmax、橫向最近距離εxmin構(gòu)成。結(jié)合助推火箭及分離殘骸彈道運(yùn)動(dòng)方程，綜合考慮分離殘骸運(yùn)動(dòng)過(guò)程中的不確定性因素影響，假設(shè)分離殘骸運(yùn)動(dòng)的某一組不確定性干擾因素為Pi=[Pi1,Pi2,…,Pin]T，殘骸落點(diǎn)距離發(fā)射點(diǎn)的縱向距離為Y(Pi)，橫向偏移距離為X(Pi)。則殘骸安全區(qū)預(yù)示模型如式(8)所示：

(8)

式中：Pr(·)表示殘骸落入某一區(qū)域的概率；P表示所有組不確定性干擾因素；Pfy1,Pfy2,Pfx1,Pfx2為對(duì)應(yīng)的目標(biāo)失效概率。殘骸安全區(qū)預(yù)示問(wèn)題的物理意義在于高效尋找安全區(qū)域的閾值εymax,εymin,εxmax,εxmin，使得殘骸落在該閾值之外的概率等于給定的目標(biāo)失效概率Pfy1,Pfy2,Pfx1,Pfx2，從而確定殘骸安全區(qū)域。

針對(duì)這一模型，常規(guī)的Monte Carlo方法需針對(duì)所有不確定干擾因素進(jìn)行大量的打靶分析，對(duì)殘骸落點(diǎn)結(jié)果進(jìn)行統(tǒng)計(jì)，確定安全區(qū)域的閾值，給出滿足目標(biāo)失效概率要求的安全區(qū)域范圍，這一方法體現(xiàn)了概率設(shè)計(jì)思想，符合實(shí)際物理意義，此外Monte Carlo方法具有非侵入性和無(wú)偏性，能夠給出精確的安全區(qū)預(yù)示范圍。然而，Monte Carlo方法需要大量調(diào)用彈道仿真模型，導(dǎo)致計(jì)算效率較低，時(shí)間成本較高，難以滿足快速設(shè)計(jì)迭代的工程需求。針對(duì)這一問(wèn)題，本文提出了一種基于優(yōu)化加點(diǎn)Kriging模型的安全區(qū)預(yù)示方法，結(jié)合Monte Carlo和代理模型的特點(diǎn)，通過(guò)序列優(yōu)化加點(diǎn)策略自適應(yīng)更新代理模型，直至滿足收斂條件，進(jìn)而結(jié)合Monte Carlo方法高效完成概率安全區(qū)分析計(jì)算。

2.2 Kriging模型

(9)

式中，f(P)=[f1(P),f2(P),…,fn(P)]T為回歸基函數(shù)，β=[β1,β2,…,βn]T為回歸系數(shù)，fT(P)β項(xiàng)對(duì)真實(shí)模型進(jìn)行全局性近似。z(P)為服從高斯分布N(0,σ2)的隨機(jī)過(guò)程，兩個(gè)訓(xùn)練樣本點(diǎn)Pi和Pj之間的協(xié)方差定義為：

Cov[z(Pi),z(Pj)]=σ2R(Pi,Pj),i,j=1,…,NT

(10)

式中，R(Pi,Pj)為相關(guān)函數(shù)，表征Pi和Pj之間的空間相關(guān)性，本文選取服從高斯分布的相關(guān)函數(shù)[11]。

利用對(duì)數(shù)似然函數(shù)，分別對(duì)式(9)的回歸系數(shù)β及式(10)中的方差σ2求導(dǎo)，給出兩者的極大似然估計(jì)：

(11)

(12)

(13)

進(jìn)而，針對(duì)某一組不確定性干擾因素Pi，Kriging模型在給出預(yù)測(cè)值的同時(shí)，還可以進(jìn)一步給出預(yù)測(cè)值的預(yù)估標(biāo)準(zhǔn)差，服從如下高斯分布：

(14)

2.3 優(yōu)化加點(diǎn)Kriging策略

為了進(jìn)一步提高安全區(qū)域分析計(jì)算的求解效率，提出了一種基于優(yōu)化加點(diǎn)Kriging模型的高效概率安全區(qū)預(yù)示方法。以求解縱向最遠(yuǎn)閾值εymax為例，根據(jù)識(shí)別的分離殘骸干擾因素分布類型及其分布參數(shù)，在整個(gè)隨機(jī)概率空間選取一定數(shù)量的Monte Carlo樣本點(diǎn)集PMC，結(jié)合縱向最遠(yuǎn)閾值的目標(biāo)失效概率Pfy1，為保證樣本點(diǎn)集PMC的相關(guān)系數(shù)不大于0.05，Monte Carlo樣本點(diǎn)數(shù)量NMC滿足如下條件：

(15)

(16)

(17)

(18)

綜上所述，分類失效概率可以表達(dá)為如下形式：

(19)

2.4 安全區(qū)預(yù)示流程

基于優(yōu)化加點(diǎn)Kriging模型的安全區(qū)預(yù)示方法流程圖如圖1所示。

圖1 基于優(yōu)化加點(diǎn)Kriging模型的安全區(qū)預(yù)示流程圖Fig.1 Flow chart of prediction method for safety control zone based on infill-sampling Kriging model

同理根據(jù)上述方法，可以確定縱向最近距離閾值εymin、橫向最遠(yuǎn)距離閾值εxmax和橫向最近距離閾值εxmin，從而給出子級(jí)殘骸概率安全區(qū)。

3 實(shí)例分析

以某型助推火箭子級(jí)殘骸安全區(qū)計(jì)算為例，對(duì)提出的基于優(yōu)化加點(diǎn)Kriging模型安全區(qū)預(yù)示方法進(jìn)行仿真驗(yàn)證。子級(jí)殘骸主要參數(shù)如表1所示，子級(jí)殘骸分離時(shí)刻及飛行過(guò)程中主要偏差量分布規(guī)律如表2所示，為簡(jiǎn)便起見(jiàn)，不考慮風(fēng)的影響。

表1 子級(jí)殘骸主要參數(shù)Tab.1 Main parameters of booster rocket′s debris

表2 子級(jí)殘骸主要偏差的分布規(guī)律Tab.2 Distribution regularity of booster rocket′s debris

3.1 與Monte Carlo方法對(duì)比分析

為了驗(yàn)證基于優(yōu)化加點(diǎn)Kriging模型的安全區(qū)預(yù)示方法的快速性和準(zhǔn)確性，首先采用Monte Carlo方法開(kāi)展助推火箭殘骸落點(diǎn)統(tǒng)計(jì)分析。取目標(biāo)失效概率為Pfy1=Pfy2=Pfx1=Pfx2=0.01，為保證Monte Carlo采樣點(diǎn)集相關(guān)系數(shù)小于0.05，根據(jù)式(15)計(jì)算出樣本點(diǎn)個(gè)數(shù)約為40 000，則選取打靶次數(shù)為40 000對(duì)助推火箭殘骸安全區(qū)進(jìn)行統(tǒng)計(jì)。同時(shí)，取相同的目標(biāo)失效概率，采用基于優(yōu)化加點(diǎn)Kriging模型的安全區(qū)預(yù)示方法進(jìn)行助推火箭殘骸落點(diǎn)統(tǒng)計(jì)，兩種方法對(duì)應(yīng)的統(tǒng)計(jì)結(jié)果對(duì)比情況見(jiàn)表3。

表3 兩種方法計(jì)算結(jié)果對(duì)比Tab.3 Prediction result of two different methods

由表3可見(jiàn)，本文提出的基于優(yōu)化加點(diǎn)Kriging模型的安全區(qū)預(yù)示方法與Monte Carlo方法相比，預(yù)示結(jié)果相對(duì)誤差小于0.2%，具有較高的計(jì)算精度，且分析次數(shù)大大降低，具有較高的計(jì)算效率。

3.2 與極限偏差疊加方法對(duì)比分析

采用傳統(tǒng)的極限偏差疊加方法，對(duì)助推火箭殘骸分離時(shí)刻的相關(guān)偏差及殘骸氣動(dòng)力系數(shù)偏差等進(jìn)行極限偏差組合，計(jì)算出來(lái)的縱向最近距離εymin為100 434.8 m，縱向最遠(yuǎn)距離εymax為271 834.848 m，橫向最近距離εxmin為243.995 m，橫向最遠(yuǎn)距離εxmax為2214.77 m。對(duì)應(yīng)的安全區(qū)預(yù)示范圍與基于優(yōu)化加點(diǎn)Kriging模型的安全區(qū)預(yù)示結(jié)果(失效概率閾值0.000 1)對(duì)比情況見(jiàn)圖2。

圖2 優(yōu)化加點(diǎn)Kriging模型方法與極限偏差疊加法對(duì)比Fig.2 Prediction result of infill-sampling Kriging model and conventional extreme deviation overlying method

從圖2可以看出，采用基于優(yōu)化加點(diǎn)Kriging模型的安全區(qū)預(yù)示方法計(jì)算的安全區(qū)范圍遠(yuǎn)遠(yuǎn)小于傳統(tǒng)極限偏差疊加方法的安全區(qū)計(jì)算結(jié)果，前者面積僅僅是后者的13%，可大大降低發(fā)射場(chǎng)安全控制實(shí)施的難度。

此外，通過(guò)改變目標(biāo)失效概率，可實(shí)現(xiàn)安全區(qū)的精細(xì)化預(yù)示，為發(fā)射場(chǎng)實(shí)施安全控制分級(jí)管理提供理論和數(shù)據(jù)支撐。圖3中給出了失效概率閾值分別為0.1、0.01和0.000 1的助推火箭殘骸安全區(qū)預(yù)示范圍，對(duì)應(yīng)的助推火箭殘骸落點(diǎn)在該區(qū)域內(nèi)的概率分別在60%、96%和99.96%以上。

圖3 不同失效概率對(duì)應(yīng)的安全區(qū)預(yù)示結(jié)果Fig.3 Safety control zone prediction result under different failure probabilities

4 結(jié)論

本文的主要研究工作和結(jié)論如下：

1)建立火箭及其分離殘骸彈道計(jì)算動(dòng)力學(xué)模型，并采用四元數(shù)方法對(duì)姿態(tài)角解算進(jìn)行處理；

2)提出基于優(yōu)化加點(diǎn)Kriging模型的安全區(qū)預(yù)示方法，結(jié)合Monte Carlo和Kriging代理模型的特點(diǎn)，給出安全區(qū)預(yù)示流程；

3)以某型助推火箭殘骸安全區(qū)計(jì)算為例進(jìn)行仿真驗(yàn)證。結(jié)果表明，基于優(yōu)化加點(diǎn)Kriging模型安全區(qū)預(yù)示方法具有較高的準(zhǔn)確性和高效性，滿足快速迭代的工程需求，相比傳統(tǒng)極限偏差疊加方法，可顯著降低安全區(qū)覆蓋面積，具有較強(qiáng)的工程應(yīng)用價(jià)值。