考慮不確定性的模塊化戰術導彈優化設計

劉鈞圣， 王剛， 王琨， 史宏博， 郭斌

(1.西北工業大學 航天學院, 陜西 西安 710072; 2.西安現代控制技術研究所, 陜西 西安 710065)

0 引言

依托于武裝直升機平臺，直升機載空地導彈具有機動性強、作戰距離遠、殺傷威力大、命中精度高等特點，是最重要、最有效的反裝甲武器之一。目前，直升機載空地導彈發展的一個方向是模塊化、系列化[1]，即彈體采用分艙段模塊化結構安裝方式，可依據軍方使用需求，針對不同制導體制的導引頭(如：圖像類、激光類、雷達類、復合制導類)與戰斗部(如：破甲殺傷型、溫壓型、侵徹型)，靈活快速地將其組裝在一起，便于部隊作戰使用[2]。模塊化導彈具有成本低、各部件功能獨立性強、便于檢測維修與故障定位等優點。

然而在模塊化導彈設計過程中，由于存在多種組合方式，使其有效載荷、質心位置均在一定范圍內浮動，受環境影響其實際飛行速度與設計值也可能不同；另一方面，在工程研制過程中，導彈彈道參數(飛行高度、彈道傾角等)設計，會根據飛行試驗后的真實氣動參數、擾動與干擾力矩等進行調整。上述不確定性因素會導致采用傳統確定性優化得到的設計結果無法滿足設計要求。國內外眾多學者曾探討和總結了考慮不確定性的大型導彈或飛行器設計方法與設計準則。陳建江等[3]建立了基于神經網路的飛航導彈穩健性優化設計模型，設計結果顯示增大助推器有助于提高不確定性因素對主級導彈射程的影響。馬英等[4]將不確定性設計理論應用于空射巡航導彈概念設計中，研究表明：設計參數波動5%，穩健性方案的射程波動幅度比確定性最優方案減小了61%，最大飛行速度波動幅度減小了57%。Jaeger等[5]針對飛機概念設計階段存在的總體參數不確定性，利用二次響應面代理模型進行了不確定性分析，結果表明：通過增大機翼展弦比與后掠角，使飛機滿足滑跑距離、航程、彈翼展長等約束的概率得到了顯著提升，但付出了稍許質量代價。目前關于小型戰術導彈的不確定性設計研究很少。劉常青等[6]以200 kg級別防空導彈為研究對象，建立了基于偏差量的導彈總體方案設計模型。基于原型彈，擬定一組設計變量并給出統計得到的概率分布，然后通過MATLAB/Simulink彈道仿真得出導彈性能及相關約束條件。仿真結果表明：導彈推力、發射初始角、助推段末速對導彈質量與射程影響較為敏感。小型戰術導彈設計時具有自身的特殊性，其有效載荷、質心位置、飛行速度等參數對其翼載荷、推重比、發動機裝藥、靜穩定度等影響較為敏感，如果拋開上述不確定性因素去設計模塊化導彈，會導致導彈實際的飛行性能偏離確定性優化最優解很多，同時設計約束也無法滿足。

針對上述問題，本文以某型模塊化直升機載空地導彈為研究對象，將穩健性優化設計理論[7-9]融入導彈總體設計原理中，建立模塊化戰術導彈穩健性優化設計模型，旨在解決滿足設計約束條件下，減小載荷質量與質心位置、飛行速度、彈道方案3類不確定性因素對模塊化導彈起飛質量影響這一問題，并依據穩健性優化結果，提出模塊化直升機載空地導彈總體參數穩健性設計準則。

1 直升機載空地導彈總體方案

某型模塊化直升機載空地導彈用于打擊地面坦克、車輛、工事及有生力量，它以裸彈的形式掛載于直升機發射架上，起飛質量不超過38 kg，射程2～10 km，需用過載≥5g，g為重力加速度。導彈總體方案如圖1所示，采用正常式X-X氣動布局、大展弦比平直矩形彈翼(發射前彈翼向后折疊，發射后彈翼按固定時序張開)，兩級雙室雙推固體火箭發動機為導彈提供助推與續航，助推發動機不與導彈分離。彈體各艙段為模塊化設計方案，導引頭艙與引戰艙可根據作戰需求進行更換。

考慮到武裝直升機對導彈起飛質量有著嚴苛的要求，本方案在滿足射程、需用過載條件下，通過減小續航發動機裝藥，使導彈在末端以無動力滑翔方式來攻擊地面目標，來減小導彈質量，速度曲線(速度v與時間t的關系)和彈道方案(發射系y軸向位移與x軸向位移的關系)分別如圖2、圖3所示。圖2中：tb為助推發動機工作結束時刻，vb為tb時刻對應的導彈速度；ts為續航發動機工作結束時刻，vs為ts時刻對應的導彈速度；te為導彈命中目標時刻，ve為te時刻對應的導彈速度，即導彈末速。圖3中：H為導彈中制導段飛行高度，θc為導彈爬升段彈道傾角，θg為導彈俯沖段彈道傾角。

圖2 速度曲線Fig.2 Velocity characteristic

圖3 彈道方案Fig.3 Flight trajectory profile

2 直升機載空地導彈穩健性優化設計

與傳統確定性優化相比，穩健性優化設計最大的特點是要進行不確定性分析，計算目標的均值與方差。不確定性分析的經典方法為蒙特卡洛仿真(MCS)，具有計算簡單、精度高等特點。應用MCS需要樣本量要足夠大，考慮到本文導彈總體設計時，迭代計算氣動力較為耗時，不可能直接將MCS應用于導彈穩健性優化設計中。為此，引入Kriging代理模型解決該問題。

構建的導彈穩健性優化設計框架如圖4所示，包括導彈總體設計、代理模型與穩健性優化三部分，整個計算程序在MATLAB環境[10]中實現。具體優化設計流程為：

1) 定義穩健性優化問題，包括設計變量、不確定性變量/參數、優化目標、約束條件等；

2) 建立導彈總體設計模型，依據定義的設計變量評估外形參數、氣動力、動力特性、彈道方程以及導彈質量。如果前后兩次導彈質量差小于0.05 kg，則認為程序收斂，輸出性能與約束條件，否則更新外形參數進行迭代；

3) 利用試驗設計生成樣本點，然后根據建立的總體設計模型對樣本點進行分析，得到一組輸入/輸出，基于該組數據選用近似方法構建代理模型，如果代理模型精度未達到指定值，則添加樣本點直至滿足精度；

4) 對代理模型進行MCS不確定性分析，計算目標函數的均值與方差，再利用多目標優化算法求解Pareto解集，根據偏好選取設計方案。

圖4 導彈穩健性優化設計框架Fig.4 Optimization framework for missile robust design

2.1 穩健性優化設計

某模塊化直升機載空地導彈標稱有效載荷mp=15 kg，標稱靜穩定度K=8.5%. 表1列出了設計變量及取值范圍，表2給出了約束定義及范圍。圖5給出了與設計變量和約束相關的參數定義，圖中：l為導彈長度，lk為質心與焦點之間的距離，lwf為彈翼前緣與尾舵前緣之間的距離，bw為彈翼展長，bf為尾舵展長，cw為彈翼弦長，AR為彈翼展弦比。

表1 設計變量的定義及取值范圍Tab.1 Definition and value range of design variables

基于前文論述的模塊化導彈不確定性因素，確定了如表3所示的3類不確定性問題。

表2 約束的定義與范圍Tab.2 Definition and range of constraints

圖5 導彈外形參數定義Fig.5 Definition of missile layout

表3 3類不確定性優化問題描述

Tab.3 Description of 3 uncertainty problems

不確定性問題涉及變量與參數概率分布類型載荷質量mp/kg服從均勻分布[13, 16]質心位置K/%服從均勻分布[6, 11]速度方案vb/(m·s-1)在設計值處,服從正態分布,標準差為8vs/(m·s-1)在設計值處,服從正態分布,標準差為8彈道方案H/m在設計值處,服從正態分布,標準差為20θc/(°)在設計值處,服從正態分布,標準差為1θg/(°)在設計值處,服從正態分布,標準差為1

該導彈穩健性優化是在滿足射程Dr=10 km、需用過載nr=5g條件下，考慮上述3類不確定性問題，通過求解表1中的設計變量，使導彈質量均值μmt與標準差σmt最小。為避免由于不確定性因素變化而導致約束無法滿足，對約束進行概率性描述，均取90%. 優化問題可描述為：

Set 表3中的3類不確定性問題，

Findx=[vb,vs,H,θc,θg,bw,bf]T，

Min [μmt,σmt]T,

采用描述性取樣策略，樣本均值與方差的計算公式為

(1)

式中：N為樣本點個數，設定為50 000；mi為第i個隨機樣本點的質量。由(1)式易知，穩健性優化設計為多目標優化問題，本文采用非支配解排序遺傳算法(NSGA-Ⅱ)[11-12]求解該問題，設定NSGA-Ⅱ優化算法的初始種群為300，種群進化代數為500.

2.2 導彈總體設計

2.2.1 外形參數

導彈外形參數包括彈長、彈徑、彈翼展長與弦長、尾舵展長與弦長。其中，彈徑參考基準導彈與戰斗部威力定為0.14 m，彈長根據基準導彈結合發動機裝藥形式、裝藥質量及裝藥密度進行估算。此處只需計算彈翼與尾舵尺寸。彈翼尺寸由(2)式確定：

(2)

式中：mm為主級(續航)導彈質量；pm為主級導彈翼載荷；S為彈翼參考面積。pm與mm將在后續推進系統與質量模塊中分別計算。尾舵面積則由需用過載nr與最大配平舵偏角(取20°)來確定。

2.2.2 氣動特性

采用渦格法(VLM)[13]計算導彈升力、焦點、誘導阻力、尾舵效率。引入普朗特- 格勞厄脫壓縮性修正因子來計算亞聲速狀態下的翼面升力。圖6(a)顯示了利用VLM建模后的導彈幾何外形，彈翼和尾舵在建模時均分為內外兩部分，內側代表與彈身的融合區，外側代表單獨的翼面，其中尾舵外側區域可以偏轉。圖6(b)顯示了該導彈在馬赫數Ma=0.3，攻角α=2°，舵偏角δ=-4°條件下的壓力云圖。

圖6 VLM對X-X布局導彈建模Fig.6 Modeling X-X configuration missile by VLM

阻力計算基于部件疊加的原理，將全彈阻力分解為彈翼阻力和彈身阻力。亞聲速彈翼阻力由誘導阻力和型阻組成，誘導阻力由VLM直接得出，型阻則通過基于XFOIL二維翼型數據的片條理論來預測[14]；彈身阻力由摩擦阻力、壓差阻力和底部阻力構成，全部根據工程估算法[15]來計算。

基于VLM的氣動力計算方法，將全彈劃分為 480個網格，利用Intel Xeon 2.80 GHz電腦運算，每計算一組氣動力(升阻力系數關于馬赫數與攻角的數據)耗時約56 s. 該方法與風洞試驗對比結果表明：當Ma在0.3～0.6，攻角在0～8°范圍時，各氣動參數最大計算誤差[16]在12%以內，可滿足概念設計階段所需。

2.2.3 推進系統

參考成熟型號，助推級發動機比沖Ib≈220 s，續航級發動機比沖Is≈190 s.

2.2.3.1 助推發動機主要參數

工作時間tb：

(3)

式中：nxmax為導彈最大軸向過載，取20g. 將助推發動機燃料質量mb與導彈起飛質量mt的比值定義為助推發動機燃料質量相對因數μb[17]：

(4)

則

(5)

式中：X為阻力；Tb為助推發動機推力；由于助推段推力遠大于阻力，可取ξ=1. 將(3)式中的tb代入(5)式，可得μb.

2.2.3.2 續航發動機主要參數

(6)

式中：v為導彈飛行速度；Cx為阻力系數；θ為彈道傾角；ρ為空氣密度；p為該時段的翼載荷。由(6)式可知，為計算Δt內推重比需知道該時間段內的翼載荷。按照可用過載大于需用過載的要求，得到p為

(7)

2.2.4 彈道方程

前面已經計算了助推發動機的燃料相對質量因數，為了獲得續航發動機燃料相對質量因數，需根據導彈射程Dr，計算所需的續航發動機燃料質量。考慮初步設計時，事先難以確定導彈質量、發動機推力和彈翼面積，故而無法求解微分方程，為此引進相對量參數μ，它表示主級導彈在t時刻所消耗的燃料相對質量因數：

(8)

(9)

式中：Cy為升力系數；x、y分別為發射系x軸、y軸方向的位移。針對本文存在無動力滑翔段，可先通過滑翔高度H、滑翔段彈道傾角θg計算出滑翔距離Dg，然后根據總射程Dr獲得有動力飛行的距離Dp，再將Dp作為(9)式積分終點進行求解，最終可得到續航發動機燃料相對質量因數μs.

2.2.5 質量評估

導彈起飛質量mt由助推器質量mb0和主級導彈質量mm組成：

mt=mb0+mm，

(10)

由于本方案助推發動機不與主級導彈不分離，因此mb0=mb. 經過推導，可得到mm的表達式為

(11)

式中：載荷質量mp包括戰斗部、彈上制導部件、電氣及直屬部件、電池等；K1和K2分別為

(12)

Ks為導彈結構相對質量因數(包括彈身與翼面結構質量因數，具體計算見參考文獻[18])，αen為發動機結構質量因數(為發動機的結構質量與燃料質量之比，αen一般比較穩定，統計值在0.6～0.7之間)。

由(4)式、(10)式和(11)式可確定導彈起飛質量：

(13)

2.3 代理模型技術

本文構建代理模型的流程為：

1) 利用拉丁超立方抽樣(LHS)[19]試驗設計方法對由表1和載荷質量及靜穩定度組成的9維度設計空間進行抽樣，生成200個初始樣本點和60個測試樣本點；

2) 生成的樣本點，利用導彈總體設計模型對樣本點進行精確分析，得到導彈質量與約束；

3) 將獲得的這一組輸入/輸出，利用Kriging模型創建代理模型。Kriging模型[20]表示為

g(u)=f(u)+z(u)，

(14)

式中：g(u)為未知Kriging模型；f(u)為已知的關于u的全局近似函數，此處選用0階模型(f(u)為常數β)；z(u)是均值為0、方差為σ2的隨機過程。與樣本點有關的協方差為

Cov(z(u(i),z(u(j)))=σ2R[R(u(i),u(j))]，

(15)

式中：R為相關矩陣；R為相關函數，本文選用高斯函數；i=1,2,…,ns;j=1,2,…,ns；ns為樣本點個數。相關矩陣R對稱正定。高斯函數表達式為

(16)

式中：θ為未知的相關參數；Nv為設計變量維度，本文Nv=7. 未知點u處的預測值(u)的表達式為

(u)=+rT(u)R-1(g-f)，

(17)

式中：g為長度為ns的列向量，為樣本點的真實響應值；f是一個長度為ns的單位列向量；rT(u)為未知點u和樣本數據之間的相關向量，其表達式為

rT(u)=[R(u,u1)R(u,u2) …R(u,uns)]T;

(18)

=(fTR-1f)-1fTR-1g.

(19)

方差估計值的表達式為

(20)

(16)式中的相關參數θ可由最大似然估計給出：

(21)

(22)

式中：N為測試樣本點個數；qsi和si分別為第i個測試樣本點所對應的精確值和響應值；σe代表了所有樣本點相對誤差在周圍的集中程度；和σe的值越接近0，表示代理模型精度越高。

5) 如果定義的上述兩個統計量均小于5%，那么代理模型構建完畢，否則返回第1步，增加50個樣本點更新Kriging模型，直至滿足收斂標準。

通過2次迭代，當樣本量達到300時，構建的代理模型可滿足擬合精度要求，如表4所示。由于計算中發現僅有ve與AR為有效約束，此處僅給出了ve與AR的代理模型。Kriging構建關于導彈質量、末速與展弦比的代理模型關于60個測試樣本點的預測結果如圖7所示。

表4 代理模型擬合精度Tab.4 Fitting precision of surrogate model

圖7 代理模型樣本點測試Fig.7 Sample points of surrogate model

3 結果與分析

為了對比，利用確定性優化算法，設計了質量最小的導彈。下面分別針對3類不確定性問題，將穩健性與確定性優化設計結果進行討論與分析。

3.1 載荷質量與質心位置不確定性

圖8給出了載荷質量與質心位置不確定時，穩健性優化得到的導彈質量Pareto前沿，同時列出了確定性優化最佳設計點P0的不確定性分析結果。由圖8可知，P0的質量標準差高達1.86 kg，說明此類導彈對質量與質心位置變化非常敏感。為了探究載荷質量與質心位置變化對導彈起飛質量與約束的影響，從 Pareto前沿中選取了3個典型的設計點P1、P2和P3進行分析，對應得到的相關參數如表5所示。通過對比，存在如下特點：

1) 從P0到P3質量方差σmt逐漸降低。根據(13)式，如果燃料質量系數越小，mp變化引起的mt變化則越小。而在射程給定的條件下，導彈滑翔距離越大，有動力飛行距離就越小，可降低燃料質量系數。因此，為了降低σmt，需要增加滑翔距離。飛行動力學表明：導彈飛行高度H越高，無動力滑翔初速vs與末速ve之比越大、阻力系數越小，滑翔距離就越遠。由表4穩健性優化結果可知，設計點P3的飛行高度最高、vs最大，從而使其σmt最低。但是，增加vs會使續航段推重比增加，導致全彈質量變大，即μmt從P0到P3逐漸增加。

圖8 載荷質量與質心位置不確定時的穩健性優化Pareto前沿Fig.8 Pareto front of RDO under payload weight and center-of-mass location uncertainty

表5 載荷質量與質心位置不確定時的穩健性優化設計結果

Tab.5 Results of RDO under payload weight and center-of-mass location uncertainty

輸出參數P0P1P2P3vb/(m·s-1)178.8182.3194.9199.9vs/(m·s-1)164.6165.5167.7170.2H/m143.2149.8160.1184.4設計θc/(°)8.927.397.317.64變量θg/(°)-15.0-15.0-15.0-15.0bw/m0.9220.9100.8780.862bf/m0.3650.3780.3810.390質量μmt/kg34.5234.6134.9535.59σmt/kg1.8611.6951.5921.532約束P{ve≥160m/s}/%33.8990.0494.0093.51P{AR≤12}/%31.32100.00100.00100.00

3) 載荷質量不確定性使得AR滿足概率的約束大幅度降低。由于載荷質量會減小，在設計時所需要的彈翼面積同步減小，使得展弦比增大，增加了展弦比進一步增大Cymax，從而提高翼載荷。根據(2)式，pm越大，彈翼面積就越小，迭代中這就使得AR越來越大，最終不滿足設計約束。根據優化結果，為了提升AR滿足約束的概率，需要降低bw. 需要注意的是，減小彈翼展長會使展弦比降低，從而降低奧斯瓦爾效率因子，增大誘導阻力，導致全彈起飛質量變大。

4) 由(5)式可知，增加助推末速vb使助推段燃料增大；由(6)式可知，增加vb可降低dv/dt，進而減小續航段推重比實現續航段燃料質量的降低。因此，從mt最小的角度來看，助推發動機與續航發動機對vb的要求是矛盾的。由表4中的Pareto最優設計結果可以看出，隨著vs的增加，一定程度增加vb可帶來更小的mt.

5) 通過上述分析，在載荷質量與質心位置不確定條件下，為增加質量穩健性與約束穩健性而采取的設計方法，必然會引起導彈起飛質量的增加。圖7所示的穩健性優化各設計點的μmt均要高于P0點，并且隨著σmt的減小，μmt逐漸增大。總得來說，ParetoP2設計方案能較好地兼顧均值和方差，與P0相比，σmt降低了19.8%，μmt增加了1.20%.

3.2 速度方案不確定性

速度方案vb和vs分別影響助推發動機質量與續航發動機推重比，它們直接關系著導彈的起飛質量。減小vs可以降低續航段推重比，從而減小續航發動機質量，但同時會降低導彈末速ve，減小滑翔距離，影響導彈質量穩健性。因此，當速度方案不確定時，需合理匹配其余總體參數，才能在性能和約束之間做出較好的折衷。

圖9 飛行速度不確定時的穩健性優化Pareto前沿Fig.9 Pareto front of RDO under velocity uncertainty

圖9給出了速度方案不確定時，穩健性優化得到的導彈質量Pareto前沿解集，同時列出了確定性優化最佳設計點P0的不確定性分析結果。從圖9可知，飛行速度不確定性引起的全彈σmt(約0.34 kg)要明顯低于載荷質量與質心位置變化引起的全彈σmt(約1.86 kg)。類似地，從Pareto前沿中選取了3個典型的點P1、P2和P3進行分析，對應得到的相關參數如表6所示。

表6 速度方案不確定時的穩健性優化設計結果Tab.6 Results of RDO under velocity uncertainty

優化結果表明：與載荷質量和質心位置不確定影響類似，為降低導彈σmt，需增加飛行高度H與滑翔段初速vs；為滿足末速約束概率，需提升vs；為滿足展弦比約束概率，需減小bw；而增加飛行高度H、vs，減小bw均會導致μmt增加。相比于確定性設計方案P0，P2的質量均值增加了1.90%，標準差降低了43%，滿足末速約束從49.79%提升至90.74%，滿足展弦比約束從58.37%提升至100%.

3.3 彈道方案不確定性

圖10 彈道方案不確定時的穩健性優化Pareto前沿Fig.10 Pareto front of RDO under trajectory uncertainty

圖10顯示了彈道方案不確定時，穩健性優化得到的質量Pareto前沿。由圖10可知，彈道方案不確定時，確定性優化設計點P0的σmt僅為0.07 kg，遠遠低于前文其他因素不確定性造成的質量波動。因此，在總體設計階段，此類導彈彈道方案在設計狀態小范圍浮動不會對其起飛質量產生顯著影響。

盡管彈道方案不確定性對導彈質量影響較小，但飛行高度與俯沖段彈道傾角均會影響無動力滑翔階段的末速，進而影響導彈翼載荷，使滿足約束的概率降低。因此在表7中，為了增加滿足末速約束的概率需要提高vs，為了滿足展弦比約束的概率需要降低bw.

表7 彈道方案不確定時的穩健性優化設計結果Tab.7 Results of RDO under trajectory uncertainty

4 結論

1) 針對模塊化導彈設計中存在的多種不確定性因素，本文以導彈起飛質量為目標，提出一種模塊化直升機載空地導彈穩健性優化設計方法。優化結果表明：相比于傳統確定性優化，穩健性優化在付出很小質量代價基礎上，可降低不確定性因素對設計目標的影響，大幅度提升滿足設計約束的概率。

2) 通過對工程研制過程中常見的3類不確定性問題分析，結果表明：載荷質量與質心位置變化對模塊化直升機載空地導彈的起飛質量影響最大，其次是速度方案變化，而彈道方案不確定性對此類導彈起飛質量幾乎沒有影響。

3) 模塊化直升機載空地導彈起飛質量與約束穩健性設計準則：增加導彈飛行高度與提升續航階段末速可降低不確定性變量對模塊化直升機載空地導彈起飛質量的影響；一定程度上增加助推段末速能降低導彈起飛質量；增加續航段末速、減小彈翼展長與增加尾舵展長提升了滿足約束的概率。

4) 本文構建的模塊化戰術導彈穩健性優化設計框架，還可用于研究彈翼、尾舵、發動機等部件模塊化產生的不確定性對導彈設計的影響，提高全構型模塊化戰術導彈設計穩健性。