基于組合優化算法的空面導彈全空域彈道優化

祝雯生，呂先敬，侯振乾，張順家，畢鵬，仲康

(1.上海機電工程研究所，上海 201109；2.空軍裝備部駐上海地區軍事代表局，上海 200030)

在現代高動態攻防對抗中，導彈以最快時間到達目標，對作戰效能有著極大的影響，如何在對飛行彈道最大高度、攻角、動壓、末速等多項約束下實現最快時間到達目標，是彈道設計的工作重點[1]。隨著載機性能的不斷提升，其飛行包線范圍越來越大，導致其掛載導彈的發射條件越來越復雜[2]。而彈道制導律中各個參數都會受到發射條件的影響，制導律需要在全空域內自適應變化，實現導彈全空域最快時間到達。因此，需要開展全空域彈道優化工作。

全空域是指根據導彈戰技指標、發射條件、目標類型、作戰場景等要素所確定的導彈飛行空域。全空域彈道優化是指在給定導彈氣動參數、推進系統參數、全彈質量特性的前提下，通過優化制導律參數，在多種約束條件下實現全空域任意發射條件下導彈飛行彈道的性能最優。由于需要對全空域下的各個發射條件進行彈道優化，優化工作量巨大。為此，需要構建全空域彈道制導律代理模型，在此基礎上開展科學、合理的彈道優化工作。代理模型技術可以通過少量的樣本點即可較為精確地預估其他樣本點的響應值，從而大幅度提高優化的效率，顯著減少計算量[3]。

國內外對彈道優化問題已有諸多研究，湯祁忠等[4]基于簡化的縱向平面彈道模型，在經典遺傳算法的基礎上，采用直接單重打靶法離散化控制變量，設計了末端多約束彈道的優化方法。李怡昕等[5]采用高斯偽譜法對變推力導彈的彈道進行優化，其仿真結果表明高斯偽譜法應用于可變推力導彈彈道優化設計的有效性。李偉喆和陳萬春[6]分別以飛行時間最短和末速最大為優化目標，采用hp自適應Radau 偽譜法求解雙脈沖中程空空導彈多段多約束的彈道優化問題，并根據求解結果分析了不同飛行條件對導彈性能的影響。韓聰聰等[7]采用基于單純形-遺傳算法，以二級彈道導彈為例，建立導彈主動段的彈道射程優化模型。優化結果表明：此復合優化算法在彈道優化應用中是可行有效的，并具有精度高、收斂速度快的特點。文獻[4-7]彈道優化工作方面取得了許多研究成果，然而對全空域范圍的彈道優化問題缺乏相關研究。

本文以空面導彈為對象，研究其在全空域范圍內的彈道優化問題，采用組合優化算法，以平均速度為優化目標，考慮高度、攻角、動壓、末速等多項約束條件，對空面導彈的飛行彈道進行優化，為全空域制導律設計提供支撐。

1 彈道優化模型

1.1 彈道模型

空面導彈的彈道模型通過求解三自由度質點彈道方程完成，根據設定的導彈和目標信息（發射條件）計算導彈的質心運動規律，其計算過程包括：彈目相對運動解算、確定標準大氣參數、制導律和平衡方程計算、氣動力和推力計算、質量質心計算、導彈運動參數和目標運動參數更新。計算時首先迭代求解平衡方程，獲得導彈飛行所需的平衡攻角及平衡舵偏角，由此計算出作用于導彈的氣動力，同時計算出當前的推力和質量，進而積分求解出導彈彈體動力學方程，最終給出隨時間變化的導彈質心運動規律[8]。

空面導彈在三維空間中攻擊目標過程中，導彈和目標的相對運動關系如圖1 所示。圖1 中Oxyz為慣性坐標系；M、T分別為導彈和目標（均視為質點）；r為導彈與目標之間的相對距離；qε、qβ分別為彈目視線高低角和方位角；vm、vt分別為導彈和目標的速度；θm、φm分別為導彈的彈道傾角和彈道偏角；θt、φt分別為目標的彈道傾角和彈道偏角。

圖1 彈目相對運動示意圖Fig.1 Schematic diagram of relative motion between missile and target

根據圖1 所示的導彈與目標之間相對運動關系，可以建立三維空間下導彈-目標的相對運動方程[9]：

導彈的動力學方程為

式中：ax、ay和az分別為導彈加速度在x、y、z這3 個方向上的分量。

導彈的運動方程[10]可由導彈速度vm在慣性坐標系Oxyz的3 個坐標軸分解得到，可表示為

式中：x˙m、y˙m、z˙m為導彈質心在慣性坐標系下的三軸速度。

同理，目標運動方程也可由目標速度vt在慣性坐標系Oxyz的3 個坐標軸分解得到。

1.2 優化控制變量

空面導彈的典型飛行過程可分為3 個階段：初制導段、中制導段和末制導段[11]，如圖2 所示。初制導段以初始高拋角射入較高的空域，形成利于中制導段的飛行態勢；中制導段繼續以高拋態勢飛行，快速進入高空，低大氣密度能夠減少氣動阻力，提高導彈射程和末速；末制導段導引頭截獲目標，導彈以指定角度（落角）攻擊目標，以保證對目標的打擊精度[12-14]。因此，初制導高拋角和落角的選擇很大程度決定了導彈的飛行彈道。為此，本文選擇高拋角θΔ和落角θf為優化控制變量。

圖2 空面導彈典型飛行過程Fig.2 Typical flight process of air-to-surface missile

1.3 性能指標

平均速度是空面導彈關鍵的性能指標。平均速度快，飛行時間短，不僅可以為載機盡快脫離戰場提供條件，還可以降低導彈慣性組件的積累誤差。本文以平均速度最大為優化目標，在全空域范圍內對空面導彈的彈道進行優化設計。

目標函數選擇為

式中：Vaverage為導彈的平均速度。

1.4 約束條件

本文主要考慮以下5 種約束條件：

1）高度約束。高空大氣密度低，如在4 5 km 高度處大氣密度僅為海平面的0.16%，因此，空面導彈通常會飛行到高空以減小氣動阻力。但如果導彈飛行高度過高，其飛行路徑將隨之增長，由于受到發動機總能量限制，速度并不會無限增加，存在上限，此時過長的飛行路徑反而會降低平均速度。因此，需要設置高度約束，即

式中：Hmax為導彈最大飛行高度。

2）攻角約束。在飛行過程中，為減少導彈能量的損耗，保證整個飛行過程狀態穩定，必須約束攻角的范圍，即

式中：αmax和αmin分別為導彈飛行過程中最大攻角和最小攻角。

3）末速約束。為了保證導彈在末制導段有足夠的能量完成對目標的有效攻擊，需要設置末速約束，即

式中：Vend為導彈的末速約束值。

4）動壓約束。為保證導彈的飛行品質和可控性并滿足結構強度要求，要求動壓在一定范圍內。因此給定動壓約束，即

式中：qmax和qmin分別為導彈飛行過程中的最大動壓和最小動壓

5）舵偏角約束。為了保證導彈飛行姿態的穩定，舵偏角的大小通常是有限的，因此給定舵偏角約束，即

式中：δmax和δmin分別為導彈的最大舵偏角和最小舵偏角。

2 彈道優化方案設計

在全空域范圍內進行彈道優化，需要將空域劃分成若干個典型點，隨后對其中每個典型點進行優化，這就導致計算量異常龐大。為此，需要構建全空域的代理模型，通過較少的樣本點較為精確地預估其他樣本點的響應值，隨后對單個樣本點進行優化，可以提升彈道優化的效率[15]。

本文建立全空域彈道2 級優化框架，如圖3 所示，其基本思想是：①將優化流程分為2 個層次，即代理模型構建和子系統彈道優化；②把全空域范圍內的導彈發射條件作為全局設計變量，全空域彈道代理模型的任務是構建導彈發射條件和優化控制變量（高拋角和落角）之間的關系；③子系統彈道優化的任務是在給定發射條件下，調整局部設計變量，在滿足約束條件下，使彈道優化的目標最優。

圖3 全空域彈道2級優化框架Fig.3 Two-level optimization framework for trajectory optimization in full airspace

2.1 代理模型

常用的代理模型有多種，如多項式響應面模型、徑向基神經網絡模型、K riging 模型和深度學習神經網絡模型等，本文采用多項式響應面模型。采用響應面模型可以構造簡單的代數表達式（多項式），考慮到空面導彈彈上計算機計算能力有限，多項式計算簡單，可滿足彈上計算機工程應用。

由于響應面模型在局部范圍內比較精確的逼近函數關系，本文將空域按照高度進行劃分，對各個高度的空域構建響應面模型，以提高代理模型的精度。代理模型的精度通常用復相關系數（R2）和均方根誤差（root meansquare error，RMSE）來衡量。其中，R2用來衡量代理模型與樣本點的近似程度，一般認為R2＆gt;0.9時，近似模型的精度可以接受；RMSE 用來衡量樣本的離散程度，一般認為小于0.2 即滿足要求[3]。

2.2 組合優化算法

目前，用于優化分析的優化算法通常可分為直接搜索算法、梯度優化算法和全局優化算法三類算法。直接搜索算法無需梯度信息，但有可能落入局部最優解；梯度優化算法優化效率高，但非常依賴初始設計點；全局優化算法搜索效率低，但能夠有效搜索全局最優解。

本文空面導彈的全空域彈道優化問題是一種具有多種約束條件的非線性優化問題，采用直接搜索算法、梯度優化算法容易陷入局部極值狀態，不能或不一定找到真正的最優解。目前的研究更傾向采用全局優化算法求解彈道優化問題，其對一些大型復雜非線性優化問題表現出良好的性能，缺點就是局部搜索效率低，最終的求解精度較差。

為此，本文提出一種基于多島遺傳算法[16]（multi-island geneticalgorithm,M IGA）和序列二次規劃法[17]（sequential quadraticprogramming,SQP）的組合優化算法，所提算法結合全局優化算法和梯度優化算法的優點，具有魯棒性強、高效、可靠、對初值不敏感等優點，可以保證全空域各個發射條件下的彈道優化都可以得到滿意的最優解。

1）M IGA 優化算法。M IGA 優化算法是一種對傳統遺傳算法的改進算法，其將傳統遺傳算法的種群進一步劃分為若干子群（島），如圖4 所示，在每個島上運用傳統遺傳算法（genetic algorithm，GA）進行子種群進化（選擇、交叉、變異等遺傳操作），并每隔一定的代數進行島間遷移，以維持種群的多樣性，避免提前收斂。M IGA 具有比傳統遺傳算法更優良的全局求解能力和計算效率，在處理大規模復雜問題時具有較大的優勢。

圖4 傳統GA和MIGA對比Fig.4 Comparison between traditional GA and MIGA

M IGA 的基本步驟為：①生成初始種群；②計算每一個個體的適應度；③按照精英保留策略選擇將要進入下一代的個體；④對種群進行選擇、交叉和變異運算，形成新的種群，并且每隔幾代挑選子群中的個體進行遷移操作；⑤新種群作為初始種群，重復第2～4 步直至滿足迭代停止條件。

2）SQP 優化算法。SQP 優化算法的核心思想是在每個迭代點x k構造一個二次規劃問題，以其解作為迭代搜索方向，并沿該方向進行一維搜索得到新的迭代點，通過反復迭代最終逼近約束優化問題的最優解。二次規劃子問題描述如下：

式中：H k為Hessian 矩陣的正定近似矩陣；hi(x k)和gj(x k)分別為等式約束和不等式約束函數；f(x k)為目標函數；d為搜索方向向量。

解式（10）得出當前迭代下搜索方向d k，從而得出新的迭代點：

式中：αk為步長參數。

隨后利用新的迭代點x k+1對H k進行更新：

式 中：s k=x k+1?x k；y k=G k+1?G k，G k為Lagrange 函數的一階導數。

3）M IGA-SQP 組 合 優 化 算 法。M IGA-SQP 組合優化算法的核心思想為：①采用M IGA 進行初步優化，得到精度較差但全局收斂性較好的次優解；②采用SQP 在搜索到的次優解附近進行局部精確搜索，以確保解的最優性。M IGA-SQP 組合優化算法兼具兩算法的優點，通過M IGA 確保可以找到一個全局性較好的優化解，避免陷入局部最優解的困境，同時又通過SQP 來彌補M IGA 的缺點，從而提高搜索速度和優化精度。

2.3 彈道優化流程

空面導彈全空域彈道優化流程如圖5 所示。

圖5 全空域彈道優化流程Fig.5 Flowchart of trajectory optimization in full airspace

1）在全空域范圍內，對發射條件（發射速度、發射高度、彈目距離等）構建樣本點，通過M IGA-SQP組合優化算法對每個樣本點進行優化。

2）根據彈道優化模型，確定初始種群和目標函數，通過M IGA 進行遺傳進化產生新一代種群，選擇其中優異的個體更新種群，經過多代遺傳得到具有較好全局收斂性的次優解。

3）以次優解作為SQP 算法的初值，植入SQP算法深度尋優，通過彈道計算得到彈道性能指標，隨后判斷結果能否達到要求的精度，若滿足則得到當前樣本點的最優解，若不滿足則重新進行SQP 優化迭代。

4）基于所有樣本點的優化結果，采用三階多項式響應面方法構造全空域彈道代理模型，建立全空域發射條件相對優化控制變量之間的擬合函數關系式。當發射條件改變后，全空域彈道代理模型給出此發射條件下的高拋角和落角。

3 算例分析

根據空面導彈全空域彈道優化流程（見圖5），進行全空域彈道優化。全空域發射條件的范圍如下：

式中：Hm0、Mam0分別為導彈發射高度和馬赫數；xt0、yt0分別為導彈發射時刻目標的x向和y向位置坐標。

為了得到具有工程意義的解，給出優化控制變量和約束條件的范圍如下：

在優化求解過程中，M IGA 算法取子種群規模選擇5，島的個數為4，總進化代數設為10，總迭代次數為200，交叉概率取0.8，變異概率取0.05，島間遷移率取0.05；SQP 算法取迭代次數為50，相對步長為0.01，收斂精度為0.001。

為了便于分析全空域范圍內彈道優化的結果，選取2 種典型發射條件的優化結果進行分析，如表1所示。為了體現所提組合算法的先進性和優勢，本文將2 種發射條件下，采用M IGA、SQP 和M IGA-SQP這3 種優化算法得到的迭代收斂情況進行對比分析，如圖6 所示。其中，矩形曲線為采用M IGA-SQP算法得到的收斂曲線，黑色矩形為M IGA 算法中每一代的最優值（共10 代），紅色矩形為SQP 算法得到的優化結果；圓點曲線為單獨采用M IGA 算法得到的收斂曲線，黑色圓點為M IGA 算法中每一代的最優值（共20 代）；三角曲線為單獨采用SQP 算法得到的收斂曲線。

表1 典型發射條件Table 1 Typical launch conditions

從圖6 中可以看出，單獨采用M IGA 算法可以得到精度較差但全局收斂性較好的解，缺點是需要多代優化才能得到滿意解（本文將代數設為20，總迭代次數為1 000）；單獨采用SQP 算法優化效率很高，但跟初值的選擇關系較大，常常陷入局部最優；采用M IGA-SQP 算法可以有效地提高收斂精度和收斂速度，具有很好的收斂性。

圖6 平均速度迭代收斂結果Fig.6 Iteration convergence results history of average velocity

如圖7 和圖8 所示2 種發射條件下優化前后彈道和速度對比曲線。從圖7 和圖8 中可以看出，優化后的高拋角和落角使得彈道的中制導段速度更快，速度分布更為合理，能量利用率更高。

圖7 2種發射條件下優化前后彈道曲線Fig.7 Trajectory before and after optimization under two launch conditions

圖8 2種發射條件下優化前后速度曲線Fig.8 Velocity curve before and after optimization under two launch conditions

表2 給出了2 種發射條件下彈道優化前后各個階段的平均速度對比。從表2 中可以發現，整個彈道優化對初制導平均速度的提升有限，因為初制導時間較短；優化主要用于提高中制導段的平均速度，此段飛行距離最遠，耗時最長，對整條彈道的影響也最大；當中制導段的平均速度提升后，給末制導段留下更多能量，末制導段平均速度隨之提高。彈道優化通過優化出合理的高拋角和落角，使得中制導段的能量和速度損失最低，從而提高整個彈道的平均速度。

表2 優化前后各個階段的平均速度對比Table 2 Comparison of average velocity for each stage before and after optim ization

在全空域范圍內構建發射條件對高拋角和落角的代理模型，并通過三階多項式函數表達，如式（15）所示。隨后，對代理模型進行誤差分析，其結果如表3 所示。其中，各個發射高度下高拋角和落角的R2基本大于0.9，RMSE 的值均較小，代理模型的精度較高，滿足工程要求。

表3 代理模型誤差分析Table 3 Error analysis for surrogate m odel

4 結 論

1）本文針對空面導彈全空域彈道優化問題，建立了多約束的彈道優化模型，并提出一種基于M IGA和SQP 的組合優化算法，采用M IGA 算法進行全局搜索優化，隨后采用SQP 算法在全局優化結果的附近進行局部尋優。算例結果表明，M IGA-SQP 組合算法有效地結合了2 種算法的優勢，克服了M IGA無法得到精確的全局最優值和SQP 對初值要求高的缺陷，算法精度相對較高，收斂速度較快，可以有效解決全空域彈道優化計算量大的問題。

2）采用多項式響應面代理模型技術，建立了全空域下的發射條件相對高拋角和落角的擬合函數關系式，為全空域制導律設計提供有力的支撐。算例結果表明，多項式響應面代理模型具有較高的精度，可以極大降低全空域彈道優化的計算量。

本文在全空域彈道優化過程中未將發動機參數（推力、工作時間等）考慮進去，在今后研究中，還需進行發動機與導彈總體的綜合優化設計，以實現導彈總體性能最優。