王 浩,張海瑞,王 堯,洪東跑,何沛昊

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

隨著武器裝備的日益發(fā)展，臨近空間高超聲速飛行器具有高速度、強(qiáng)突防、強(qiáng)生存等優(yōu)勢，是世界主要軍事大國的研究熱點(diǎn)之一。與常規(guī)彈道導(dǎo)彈相比，新型高超聲速飛行器具有大空域，寬速域的特點(diǎn)，其飛行環(huán)境更為惡劣、約束條件更為嚴(yán)格、不確定性影響更為顯著。通用航空飛行器(Common Aero Vehicle， CAV)是一種機(jī)動靈活性極高的臨近空間高超聲速飛行器，具有典型的代表意義，CAV飛行器既可以從空間平臺投放，也可以利用固體火箭發(fā)動機(jī)作為助推級從地面發(fā)射，具有良好的通用性和適應(yīng)性[1]。

常規(guī)射程評估方法通過偏差組合考慮最大射程的最小值，以此包絡(luò)預(yù)定射程[2]，由于忽略了系統(tǒng)固有的可變性以及模型的不確定性，在新型高超聲速飛行器射程評估中具有一定的局限性和主觀性，評估精度較差，同時難以實(shí)現(xiàn)射程裕度的量化分析。

針對新型高超聲速飛行器的特點(diǎn)，綜合考慮不確定性因素影響，開展飛行器射程概率評估具有重要意義，不僅可以作為衡量射程能力的重要指標(biāo)，還可以為射程裕度分析提供定量依據(jù)。

在工程應(yīng)用中，通常采用Monte Carlo方法對射程進(jìn)行評估。Monte Carlo方法具有非侵入性和無偏性[3]，適用于求解隱式非線性問題，在工程中有著廣泛的應(yīng)用。然而，Monte Carlo方法需要調(diào)用大量的仿真計(jì)算模型，從而消耗巨大的時間成本和計(jì)算資源。為了進(jìn)一步提升計(jì)算效率，近年來國內(nèi)外學(xué)者結(jié)合代理模型和Monte Carlo方法開展了較多的研究工作。文獻(xiàn)[5]結(jié)合Kriging模型的特點(diǎn)提出了主動學(xué)習(xí)Kriging方法，通過U函數(shù)序列加點(diǎn)，具有很高的計(jì)算效率。文獻(xiàn)[6]進(jìn)一步將主動學(xué)習(xí)Kriging和子集模擬方法相結(jié)合，提出了AK-SS方法處理小失效概率隱式非線性問題。

本文以某CAV構(gòu)型為研究對象，利用固體火箭發(fā)動機(jī)助推發(fā)射，建立了氣動、推進(jìn)、質(zhì)量、彈道學(xué)科分析模型，綜合考慮不確定性因素影響，建立飛行器射程概率評估模型，通過引入主動學(xué)習(xí)方法和Kriging代理模型，對CAV構(gòu)型飛行器射程進(jìn)行概率評估。

1 飛行器射程評估建模

1.1 飛行器總體基準(zhǔn)方案

以固體火箭發(fā)動機(jī)助推CAV構(gòu)型飛行器方案為基準(zhǔn)，假設(shè)飛行器從發(fā)射陣地點(diǎn)火起飛，經(jīng)火箭發(fā)動機(jī)助推至一定高度后頭體分離，進(jìn)而，飛行器開始機(jī)動滑翔，依靠氣動力控制實(shí)現(xiàn)跳躍滑翔，最后，在距離目標(biāo)較近的時刻，實(shí)現(xiàn)彈道下壓，進(jìn)入末制導(dǎo)段，俯沖直至命中目標(biāo)，確保飛行器的精確打擊。任務(wù)剖面可劃分為助推段、滑翔段以及末制導(dǎo)段,飛行器的任務(wù)剖面如圖1所示[7]。

圖1 新型臨近空間飛行器任務(wù)剖面

由于末制導(dǎo)階段射程與其他兩段相比較小，只需要考慮助推段和滑翔段射程即可。此外，高超聲速飛行器的射程與設(shè)計(jì)方案的氣動、推進(jìn)、質(zhì)量、彈道4個學(xué)科密切相關(guān)，給出氣動、推進(jìn)、質(zhì)量、彈道4個子學(xué)科與飛行器射程關(guān)系，射程仿真計(jì)算流程方塊圖如圖2所示[8]。

1.2 學(xué)科分析模型

氣動學(xué)科分析模型是根據(jù)飛行器總體基準(zhǔn)方案的外形布局、氣動學(xué)科設(shè)計(jì)變量以及飛行狀態(tài)參數(shù)等計(jì)算氣動力系數(shù)，包括阻力系數(shù)Cd和升力系數(shù)Cl。針對CAV構(gòu)型滑翔彈頭方案，采用文獻(xiàn)[9]氣動分析數(shù)據(jù)，結(jié)合總體方案滑翔高度及速度特點(diǎn)，選取典型滑翔高度為30 km，建立關(guān)于攻角α和飛行馬赫數(shù)Ma的阻力系數(shù)和升力系數(shù)矩陣[10]，之后，通過線性插值計(jì)算某一攻角和馬赫數(shù)條件下的氣動力系數(shù)。

圖2 射程仿真計(jì)算流程

推進(jìn)學(xué)科分析模型是根據(jù)給定的發(fā)動機(jī)裝藥形式和裝藥參數(shù)，進(jìn)行熱力計(jì)算、內(nèi)彈道性能計(jì)算、質(zhì)量計(jì)算等，最終給出壓強(qiáng)-時間曲線p(t)、推力-時間曲線T(t)以及秒耗量-時間曲線m(t)等?；谌紵覂?nèi)壓力均一假設(shè)，推進(jìn)學(xué)科分析采用零維內(nèi)彈道計(jì)算方程組進(jìn)行計(jì)算，具有計(jì)算速度快，計(jì)算精度適中的特點(diǎn)[11]。推進(jìn)方案采用單室雙推力固體火箭發(fā)動機(jī)，其藥柱包括助推段和續(xù)航段兩部分，從而在一定程度上彌補(bǔ)固體火箭發(fā)動機(jī)推力不可調(diào)節(jié)以及工作時間短的弱點(diǎn)。

質(zhì)量學(xué)科分析模型主要包括CAV構(gòu)型飛行器質(zhì)量、發(fā)動機(jī)推進(jìn)劑質(zhì)量、發(fā)動機(jī)結(jié)構(gòu)質(zhì)量、以及尾翼、控制機(jī)構(gòu)等結(jié)構(gòu)質(zhì)量。發(fā)動機(jī)推進(jìn)劑質(zhì)量可通過推進(jìn)學(xué)科質(zhì)量計(jì)算給出，發(fā)動機(jī)結(jié)構(gòu)質(zhì)量按照展開型質(zhì)量方程計(jì)算[2]，包括筒體質(zhì)量、前后封頭質(zhì)量、絕熱層和包覆層質(zhì)量、噴管質(zhì)量以及點(diǎn)火器質(zhì)量等結(jié)構(gòu)質(zhì)量。尾翼、控制機(jī)構(gòu)等質(zhì)量，在方案設(shè)計(jì)階段可按相關(guān)成熟型號確定且保持不變。在發(fā)動機(jī)工作的助推階段，飛行器質(zhì)量隨發(fā)動機(jī)推進(jìn)劑的消耗逐漸減小；在無動力的滑翔階段，忽略防熱涂層的燒蝕作用，飛行器質(zhì)量恒為CAV構(gòu)型質(zhì)量。

彈道學(xué)科分析模型是在其他學(xué)科分析結(jié)果的基礎(chǔ)上，根據(jù)彈道運(yùn)動方程實(shí)現(xiàn)飛行器自發(fā)射至末制導(dǎo)交接班的全過程飛行仿真。本文將彈道學(xué)科模型處理為軌跡優(yōu)化問題，以飛行器射程最大為優(yōu)化目標(biāo)，在初始約束、終端約束以及控制約束等因素影響下，利用hp自適應(yīng)偽譜法實(shí)現(xiàn)彈道優(yōu)化[12]，在軌跡優(yōu)化過程中，氣動分析模型通過插值計(jì)算直接嵌入至彈道方程中。hp自適應(yīng)偽譜法融合了Radau偽譜法與hp型有限元法的優(yōu)勢，在若干離散點(diǎn)上，通過多項(xiàng)式近似替代微分方程約束，從而將最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題，具有很高的計(jì)算精度，其優(yōu)勢在于當(dāng)某一離散區(qū)間的計(jì)算精度不滿足要求時，可以自適應(yīng)調(diào)整該區(qū)間內(nèi)的配點(diǎn)數(shù)h和多項(xiàng)式階次p，提高多項(xiàng)式近似的精度和效率。hp自適應(yīng)偽譜法的計(jì)算流程如圖3所示。

1.3 飛行器射程概率評估模型

針對新型高超聲速飛行器大空域，寬速域的特點(diǎn)，在學(xué)科分析模型基礎(chǔ)上，綜合考慮氣動、推進(jìn)、質(zhì)量學(xué)科的不確定性因素影響，將不確定性因素注入飛行器射程仿真模型。

假設(shè)飛行器飛行過程中的不確定性因素為X=[x1,x2,…,xn]T，給定目標(biāo)射程為d0，則在不確定性干擾因素影響下，飛行器最大射程滿足給定目標(biāo)射程d0定義為可信射程D，即

D={d|d(X)-d0≥0}

(1)

式(1)中，d為某一不確定性因素條件下的最大射程。在不確定性因素影響下，飛行器最大射程在某一區(qū)間內(nèi)呈現(xiàn)一定的概率分布，進(jìn)而，可信射程出現(xiàn)的概率定義為射程的可信程度R，即射程概率評估模型為

R=Pr(d(X)-d0≥0)

(2)

針對這一模型，通過工程經(jīng)驗(yàn)和相關(guān)成熟型號數(shù)據(jù)確定不確定性干擾因素的分布類型及其分布參數(shù)，進(jìn)而通過Monte Carlo方法開展新型高超聲速飛行器射程概率評估。然而在工程應(yīng)用中，Monte Carlo方法需要調(diào)用大量的射程仿真模型，其中彈道學(xué)科模型為軌跡優(yōu)化問題，導(dǎo)致射程概率評估效率較低。針對這一問題，提出了一種基于主動學(xué)習(xí)Kriging的飛行器射程概率評估方法，利用主動學(xué)習(xí)和Kriging模型特點(diǎn)，通過序列加點(diǎn)策略自適應(yīng)更新Kriging模型，進(jìn)而結(jié)合Monte Carlo方法高效開展射程概率評估。

圖3 hp自適應(yīng)偽譜法的計(jì)算流程

2 主動學(xué)習(xí)Kriging模型

2.1 Kriging模型

(3)

式(3)中，f(X)=[f1(X),f2(X),…,fn(X)]T為回歸基函數(shù)，n為回歸基函數(shù)的個數(shù)，β=[β1,β2,…,βn]T為對應(yīng)的回歸系數(shù)，該項(xiàng)對模型進(jìn)行全局性近似；z(X)為服從高斯分布N(0,σ2)的隨機(jī)過程，該項(xiàng)對模型進(jìn)行一定程度的修正[15]。

利用對數(shù)似然函數(shù)，給出回歸系數(shù)β及方差σ2的極大似然估計(jì)

(4)

(5)

其中，F(xiàn)=[f(X1),f(X2),…,f(XNT)]T是f(X)的向量，R為相關(guān)矩陣[16]。

在某一未知不確定性因素影響下，Kriging模型在給出預(yù)測射程的基礎(chǔ)上，還可以進(jìn)一步給出預(yù)測射程的預(yù)估標(biāo)準(zhǔn)差，服從如下高斯分布

(6)

式(6)中，μD(X)為Kriging的預(yù)測射程，σD(X)為預(yù)測射程的預(yù)估標(biāo)準(zhǔn)差。

2.2 主動學(xué)習(xí)策略

根據(jù)識別的不確定性干擾因素分布類型及其分布參數(shù)，在整個不確定性空間選取一定數(shù)量的Monte Carlo樣本點(diǎn)集XMC，樣本點(diǎn)數(shù)量Nc要求滿足相關(guān)系數(shù)小于0.05，即

(7)

根據(jù)Kriging模型服從高斯分布的特性，樣本點(diǎn)集XMC中任意點(diǎn)分類失效概率可表達(dá)為如下：

(8)

其中，U函數(shù)為學(xué)習(xí)函數(shù)，當(dāng)U函數(shù)較小時，分類失效概率π較大，在當(dāng)前Kriging模型下，該樣本點(diǎn)具有較大不確定性。進(jìn)而，在樣本點(diǎn)集XMC中選取U函數(shù)最小值對應(yīng)的樣本點(diǎn)作為序貫新增訓(xùn)練點(diǎn)，即選取最大分類失效概率的樣本點(diǎn)。

通過上述主動學(xué)習(xí)策略，序列加點(diǎn)，自適應(yīng)更新Kriging模型，若在迭代過程中，U函數(shù)的最小值大于2時，此時樣本點(diǎn)集XMC的最大分類失效概率小于Φ(-2)≈0.022 8，滿足射程概率評估精度要求，判定收斂?；谥鲃訉W(xué)習(xí)Kriging模型的飛行器射程概率評估計(jì)算步驟如下：

1) 根據(jù)識別的飛行過程不確定性干擾因素分布類型及分布參數(shù)產(chǎn)生樣本點(diǎn)集XMC。

4) 利用Kriging模型計(jì)算U(XMC)，選取最小值min(U(XMC))對應(yīng)的樣本點(diǎn)為新增訓(xùn)練點(diǎn)Xnew。

5) 若min(U(XMC))>2，Kriging模型滿足精度要求，完成飛行器射程概率評估計(jì)算；否則，繼續(xù)步驟6。

3 數(shù)值仿真

以固體火箭發(fā)動機(jī)助推CAV構(gòu)型飛行器方案為基準(zhǔn)，CAV構(gòu)型飛行器質(zhì)量為907.18 kg，氣動面積為0.483 87 m2，目標(biāo)射程為1 300 km。若在不確定性因素影響下射程可信程度大于95%，即滿足設(shè)計(jì)要求；否則，修改論證方案直至滿足射程概率評估條件。主要參數(shù)如表1所示。

表1 基線方案主要標(biāo)稱參數(shù)

忽略不確定性因素影響，利用零維內(nèi)彈道方程組給出標(biāo)稱狀態(tài)下壓力-時間曲線p(t)以及推力-時間曲線T(t)，如圖4所示。進(jìn)而，給定彈道終端約束為終端高度20 km，終端速度1 000 m/s，終端彈道傾角-23°，利用hp自適應(yīng)偽譜法給出標(biāo)稱狀態(tài)下射程最大的縱程-高程曲線d(h)以及速度-時間曲線v(t)，如圖5所示。

圖4 發(fā)動機(jī)壓力-時間p(t)和推力-時間T(t)曲線

通過仿真計(jì)算，標(biāo)稱狀態(tài)下最大射程達(dá)到了1 510 km，與目標(biāo)射程1 300 km相比預(yù)留210 km裕度，為驗(yàn)證上述射程裕度的合理性，開展射程裕度精細(xì)化評估。結(jié)合飛行器彈道特點(diǎn)以及工程經(jīng)驗(yàn)，考慮如下不確定性因素，包括助推段發(fā)動機(jī)推力偏斜θ、助推段發(fā)動機(jī)推力系數(shù)偏差ΔP，氣動升力系數(shù)偏差ΔCl，氣動阻力系數(shù)偏差ΔCd以及滑翔彈頭質(zhì)量偏差Δm，其不確定性干擾因素分布類型及分布參數(shù)值如表2所示。

圖5 彈道的縱程-高程d(h)以及速度時間v(t)曲線

不確定性因素分布類型分布參數(shù)θ/(″)正態(tài)分布N(0,1.5)ΔP均勻分布U(-0.03,0.03)ΔCl均勻分布U(-0.15,0.15)ΔCd均勻分布U(-0.15,0.15)Δm/kg正態(tài)分布N(0,3)

根據(jù)識別的不確定性因素模型隨機(jī)選取10 000個Monte Carlo樣本點(diǎn)，利用拉丁超立方采樣方法抽取21個訓(xùn)練點(diǎn)構(gòu)建初始Kriging模型，進(jìn)而采用主動學(xué)習(xí)Kriging方法開展射程概率評估，通過154次序列加點(diǎn)，自適應(yīng)更新代理模型，滿足收斂條件。進(jìn)一步給出加點(diǎn)射程響應(yīng)與目標(biāo)射程之差的迭代歷程如圖6所示。

通過圖6表明：基于主動學(xué)習(xí)Kriging的飛行器射程概率評估方法趨于選擇目標(biāo)射程附近的樣本點(diǎn)以及不確定性因素較大的樣本點(diǎn)，因而，Kriging模型在目標(biāo)射程附近具有較高的精度，滿足射程概率評估需求，為進(jìn)一步驗(yàn)證方法的準(zhǔn)確性和有效性，采用常規(guī)的Monte Carlo方法開展射程概率評估，評估結(jié)果如表3所示，其中計(jì)算仿真平臺的處理器為Intel(R) Core(TM) i7-6700U CPU @3.40GHz,內(nèi)存4.00GB。

圖6 加點(diǎn)射程響應(yīng)與目標(biāo)射程之差的歷程

由表3可見，基于主動學(xué)習(xí)Kriging的飛行器射程可信程度與Monte Carlo方法結(jié)果的相對誤差不足0.5%，具有很高的計(jì)算精度。同時，本文方法調(diào)用射程仿真次數(shù)遠(yuǎn)小于Monte Carlo方法的分析次數(shù)，具有較高的計(jì)算效率，滿足方案階段快速設(shè)計(jì)迭代的工程需求。此外，分析結(jié)果表明，飛行器基線方案的射程可信程度大于95%，滿足射程指標(biāo)需求，實(shí)現(xiàn)了射程裕度的量化分析計(jì)算。

表3 飛行器射程概率評估結(jié)果

4 結(jié)論

針對臨近空間高超聲速飛行器大空域、寬速域的特點(diǎn)，在建立氣動、推進(jìn)、質(zhì)量、彈道學(xué)科分析模型的基礎(chǔ)上，通過考慮飛行過程不確定性因素的影響，建立了飛行器射程概率評估模型，實(shí)現(xiàn)了射程裕度的定量評估。在此基礎(chǔ)上，通過引入主動學(xué)習(xí)策略和Kriging代理模型提升了射程評估效率，結(jié)合數(shù)值仿真結(jié)果驗(yàn)證了該方法的準(zhǔn)確性和有效性，減少了射程評估的計(jì)算時間，為飛行器總體方案射程評估提供定量依據(jù)，具有較強(qiáng)的工程實(shí)踐意義。