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

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

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

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

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

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

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

1 飛行器射程評估建模

1.1 飛行器總體基準方案

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

圖1 新型臨近空間飛行器任務剖面

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

1.2 學科分析模型

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

圖2 射程仿真計算流程

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

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

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

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

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

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

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

(1)

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

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

(2)

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

圖3 hp自適應偽譜法的計算流程

2 主動學習Kriging模型

2.1 Kriging模型

(3)

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

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

(4)

(5)

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

在某一未知不確定性因素影響下，Kriging模型在給出預測射程的基礎上，還可以進一步給出預測射程的預估標準差，服從如下高斯分布

(6)

式(6)中，μD(X)為Kriging的預測射程，σD(X)為預測射程的預估標準差。

2.2 主動學習策略

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

(7)

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

(8)

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

通過上述主動學習策略，序列加點，自適應更新Kriging模型，若在迭代過程中，U函數(shù)的最小值大于2時，此時樣本點集XMC的最大分類失效概率小于Φ(-2)≈0.022 8，滿足射程概率評估精度要求，判定收斂。基于主動學習Kriging模型的飛行器射程概率評估計算步驟如下：

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

4) 利用Kriging模型計算U(XMC)，選取最小值min(U(XMC))對應的樣本點為新增訓練點Xnew。

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

3 數(shù)值仿真

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

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

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

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

通過仿真計算，標稱狀態(tài)下最大射程達到了1 510 km，與目標射程1 300 km相比預留210 km裕度，為驗證上述射程裕度的合理性，開展射程裕度精細化評估。結合飛行器彈道特點以及工程經(jīng)驗，考慮如下不確定性因素，包括助推段發(fā)動機推力偏斜θ、助推段發(fā)動機推力系數(shù)偏差ΔP，氣動升力系數(shù)偏差ΔCl，氣動阻力系數(shù)偏差ΔCd以及滑翔彈頭質量偏差Δ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ù)識別的不確定性因素模型隨機選取10 000個Monte Carlo樣本點，利用拉丁超立方采樣方法抽取21個訓練點構建初始Kriging模型，進而采用主動學習Kriging方法開展射程概率評估，通過154次序列加點，自適應更新代理模型，滿足收斂條件。進一步給出加點射程響應與目標射程之差的迭代歷程如圖6所示。

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

圖6 加點射程響應與目標射程之差的歷程

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

表3 飛行器射程概率評估結果

4 結論

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