劉 欣，楊 濤，張青斌

（國防科學(xué)技術(shù)大學(xué) 航天與材料工程學(xué)院，長沙410073）

助推－滑翔導(dǎo)彈通過火箭助推加速，依靠滑翔彈頭的高升阻比外形，利用空氣動力控制飛行軌跡，實現(xiàn)遠距離的非彈道式再入機動飛行.與彈道導(dǎo)彈相比，由于其在增大射程、突防、再入段機動等方面具有優(yōu)勢，助推－滑翔導(dǎo)彈逐漸成為近來的研究熱點[1].

彈道優(yōu)化設(shè)計是助推－滑翔導(dǎo)彈總體設(shè)計的關(guān)鍵技術(shù)之一.助推－滑翔導(dǎo)彈的彈道可初步分為主動段、自由飛行段、再入段，目前國內(nèi)外相關(guān)研究集中在滑翔彈頭再入彈道優(yōu)化設(shè)計方面，針對導(dǎo)彈全彈道優(yōu)化的研究較少.文獻[2]利用近似方法對全球到達高超聲速飛行器的彈道特性進行了分析，文獻[3]利用序列二次規(guī)劃法對助推－滑翔導(dǎo)彈縱向全彈道進行最大射程優(yōu)化.由于助推－滑翔導(dǎo)彈彈道的復(fù)雜性，利用解析法進行總體參數(shù)分析十分困難，且誤差很大，結(jié)合各種優(yōu)化技術(shù)對助推－滑翔導(dǎo)彈的彈道特性與總體參數(shù)進行分析是一種趨勢.

本文利用Gauss偽譜法對全彈道進行優(yōu)化設(shè)計，并分析了起飛質(zhì)量、主動段終端傾角等參數(shù)對導(dǎo)彈射程的影響.

1 助推器參數(shù)估算

為了完成主動段彈道的計算，建立起飛質(zhì)量與主動段終點參數(shù)乃至導(dǎo)彈射程的聯(lián)系，需要對助推器的基本參數(shù)進行估算.

假設(shè)助推－滑翔導(dǎo)彈采用單級固體火箭助推.若初始速度為零，根據(jù)齊奧爾柯夫斯基公式，單級火箭的理想速度增量為[4]

式中：u′e為有效排氣速度，數(shù)值上等于真空比沖Ip；mf為載荷質(zhì)量，m0為火箭起飛質(zhì)量，m0－mf為助推器質(zhì)量；結(jié)構(gòu)系數(shù)ε為助推器結(jié)構(gòu)質(zhì)量與助推器質(zhì)量之比.考慮飛行過程中由引力、氣動力和大氣靜壓引起的飛行器速度損失Δv，則要使火箭主動段終點速度達到vf，火箭的理想速度增量應(yīng)為vid＝vf＋Δv，可求得起飛質(zhì)量與目標(biāo)速度、載荷質(zhì)量的關(guān)系：

依據(jù)文獻[5]中的相關(guān)論述，可取u′e＝Ip＝2 800m/s，ε＝0.1，Δv＝0.2vf.對于給定的彈頭，式（2）建立了起飛質(zhì)量與主動段終點速度的關(guān)系.

假設(shè)助推系統(tǒng)的任務(wù)是將質(zhì)量為907.8kg的滑翔彈頭[6]助推至速度4 500m/s.根據(jù)式（2）可求得起飛質(zhì)量m0＝18 013kg，進一步可計算出助推器裝藥質(zhì)量為mp＝（m0－mf）（1－ε）＝15 394.7kg.通常導(dǎo)彈設(shè)計時需對主動段重推比（ν0＝m0g/（mpIp/t1））加以選擇.參考一般固體火箭取值，重推比可取為ν0＝0.5，據(jù)此可求出助推器工作時間t1＝s，燃料質(zhì)量流量平均推力Fp＝Ip＝353.08kN.

2 助推－滑翔導(dǎo)彈彈道優(yōu)化模型

彈道優(yōu)化設(shè)計的目的是獲取一條或一類滿足約束和性能指標(biāo)條件的最優(yōu)彈道.彈道優(yōu)化問題可描述為一般的最優(yōu)控制問題，即在時間區(qū)間[t0，tf]中，尋找最優(yōu)控制變量u（t），最小化性能指標(biāo)J，并使?fàn)顟B(tài)變量x（t）、初始時間t0、終端時間tf滿足運動微分方程約束以及端點約束和過程約束.下面結(jié)合助推－滑翔導(dǎo)彈的特點，建立運動微分方程，確定約束條件及性能指標(biāo).

1）性能指標(biāo).

助推－滑翔導(dǎo)彈彈道復(fù)雜，宜分段進行彈道優(yōu)化.滑翔彈頭的主要優(yōu)勢即在于增大射程，總體設(shè)計時比較關(guān)心其射程覆蓋范圍，以Lf表示彈道終點處射程，可將滑翔彈頭飛行段的性能指標(biāo)取為J＝min（－Lf），即使射程最大.

主動段的主要作用是加速滑翔彈頭，其射程在全彈道射程中所占比例很小，可將性能指標(biāo)取為J＝－min[（gH＋v2/2）｜tf]，即主動段終點機械能最大，以利于滑翔彈頭飛行段增大射程.

2）狀態(tài)方程.

總體性能分析時主要關(guān)注導(dǎo)彈縱平面的運動，忽略地球自轉(zhuǎn)和非球形影響，建立助推－滑翔導(dǎo)彈縱平面運動方程：

式中：H為飛行器距離地面的高度；L為射程；Θ為當(dāng)?shù)厮俣葍A角；v為飛行器飛行速率；m為飛行器質(zhì)量；e為發(fā)動機燃料質(zhì)量流量；FX，F(xiàn)Y分別為氣動阻力、升力；r為飛行器到地心的距離；μ為地球引力系數(shù)；φp為推力與速度方向的夾角.

在優(yōu)化問題中取狀態(tài)變量x＝（HLvΘ m）T，主動段通過改變推力方向控制彈道，取控制量為φp；滑翔彈頭段通過氣動舵來改變彈道，而在僅考慮縱平面運動的情況下，飛行器動力取決于飛行器的速度、高度及攻角α，則取控制量為攻角α.

3）約束條件.

彈道設(shè)計時主要考慮飛行過程約束、終端參數(shù)約束及控制量約束.飛行過程約束，也就是飛行器在飛行過程中彈道參數(shù)必須滿足的約束.為確保飛行器安全飛行，通常要考慮熱流、法向過載、動壓等因素對彈道的約束.具體有，法向過載約束：n＝，動壓約束駐點熱流約束＝ksρ0.5v3.15＜max.其中，ny，max，qmax，max分別為飛行中所允許的法向過載、動壓和熱流的最大值.

終端參數(shù)約束是指飛行器在彈道端點需滿足的條件，即最優(yōu)控制問題中的邊界條件.主動段彈道起點參數(shù)取為導(dǎo)彈垂直上升段終點狀態(tài)，滑翔彈頭飛行段起點參數(shù)應(yīng)當(dāng)與主動段終點參數(shù)一致，滑翔彈頭飛行段末端通常有落地速度與落地彈道傾角的要求，可?。ǎ╰f）＝－60°，v（tf）＝1 000m/s.

另外，攻角、推力與速度方向的夾角等控制量的幅值和變化率不能超過限制值，可?。粒?5°，｜φp｜＜5°，｜｜＜1 （°）/s，｜p｜＜1 （°）/s.

3 Gauss偽譜法優(yōu)化基本原理

助推－滑翔導(dǎo)彈彈道形式復(fù)雜，飛行約束多，其彈道優(yōu)化是一個難點.Gauss偽譜法（Gauss Pseudospectral Method，GPM）是一種適合解決多階段多約束最優(yōu)控制問題的數(shù)值優(yōu)化方法，它將狀態(tài)變量和控制變量在一系列Legendre Gauss（LG）點上離散，構(gòu)造全局插值多項式來近似狀態(tài)變量和控制變量，經(jīng)過一系列處理，最終將最優(yōu)控制問題轉(zhuǎn)化為具有一系列代數(shù)約束的參數(shù)優(yōu)化問題，即非線性規(guī)劃問題（NLP），并進行數(shù)值求解.GPM離散最優(yōu)控制問題的主要流程如下[7，8].

1）時域變換.

由于Gauss偽譜法的配點都分布在區(qū)間[－1，1]，因此求解中一般對時間區(qū)間進行變換，通過引入變換：

將最優(yōu)控制問題的時間區(qū)間由t∈[t0，tf]轉(zhuǎn)換至τ∈[－1，1]，并將原問題按時間歸一化.

2）全局插值多項式近似狀態(tài)與控制變量.

Gauss偽譜法的離散點κ＝｛τ1，…，τK｝為K階Legendre－Gauss點，即K階Legendre多項式的根PK（τ），其中，

Legendre－Gauss點分布在區(qū)間（－1，1），再增加一個點，τ0＝－1，作為區(qū)間[－1，1）上的K＋1個插值點，構(gòu)成K＋1個Lagrange插值多項式li（τ）（i＝0，…，K），并以此為基函數(shù)近似狀態(tài)變量和控制變量：

3）動力學(xué)方程約束轉(zhuǎn)換為代數(shù)約束.

偽譜方法中，狀態(tài)變量由全局插值多項式近似，狀態(tài)變量的導(dǎo)數(shù)可通過求導(dǎo)來近似，從而將動力學(xué)微分方程約束轉(zhuǎn)換為代數(shù)約束，即

式中，τk為LG點.微分矩陣D∈RK×（K＋1）可離線確定，將式（8）代入動力學(xué)方程：

于是，將最優(yōu)控制問題的動力學(xué)微分方程約束轉(zhuǎn)換為代數(shù)約束：

4）離散條件下的性能指標(biāo).

將Bolza型性能指標(biāo)函數(shù)J＝Φ（X0，t0，Xf，tf）＋中的積分項用Gauss積分來近似，得到Gauss偽譜方法中的性能指標(biāo)函數(shù)：

5）離散條件下的終端狀態(tài)約束.

最優(yōu)控制問題往往包含終端狀態(tài)約束，根據(jù)動力學(xué)方程有：

將終端狀態(tài)約束條件離散并用Gauss積分來近似，可得：

6）多階段不連續(xù)問題的處理.

對于助推－滑翔導(dǎo)彈彈道優(yōu)化這樣的分段不連續(xù)最優(yōu)控制問題，提出了多范圍偽譜技術(shù)，即將問題分成若干子區(qū)間，在每個小區(qū)間上分別進行狀態(tài)變量和控制變量的離散，在分段點處設(shè)置連接點，同時加上連接點約束條件.例如在某一點τ1∈[τ0，τ2]不連續(xù)，記

則連接點約束條件為

基于上述的數(shù)值近似方法，原連續(xù)最優(yōu)控制問題被離散，并轉(zhuǎn)換為非線性規(guī)劃（NLP）問題.轉(zhuǎn)換所得的NLP問題可采用數(shù)值方法求解，其中序列二次規(guī)劃（SQP）是一種較好的選擇，SQP算法較為成熟，在此不做介紹.BENSON和 HUNTINGTON等人的研究[9，10]證明了經(jīng)GPM參數(shù)化得到的NLP問題的 KKT（Karush－Kuhn－Tucker）條件與離散的Hamiltonian邊值問題（HBVP）最優(yōu)性條件之間的等價性.除此之外，GPM的主要優(yōu)點是算法以指數(shù)速率收斂，收斂速率快，且能在離散點較少時得到較好的精度.

4 優(yōu)化算例

1）最大射程彈道優(yōu)化.

助推－滑翔彈道的全彈道可劃分為主動段和滑翔彈頭飛行段，導(dǎo)彈經(jīng)主動段彈道飛行，于關(guān)機點處滑翔彈頭與助推級分離，分離后的彈頭經(jīng)自由飛行段飛行，再入大氣層后依靠彈頭的高升阻比氣動性能進行滑翔機動飛行，最終以某落速和落角攻擊目標(biāo).利用GPM對多階段不連續(xù)問題的處理辦法，主動段以關(guān)機點機械能最大為性能指標(biāo)，以主動段關(guān)機點參數(shù)作為滑翔彈頭飛行段起點參數(shù)，滑翔彈頭段以射程最大為性能指標(biāo)，分段進行彈道優(yōu)化.

由于滑翔彈頭飛行段射程占總射程的絕大部分，全彈道射程主要取決于主動段終點的運動參數(shù)及滑翔彈頭的氣動特性.在助推器與滑翔彈頭參數(shù)確定的情況下，主動段終端參數(shù)由主動段控制規(guī)律決定.假設(shè)Hf，vf，Θf分別為主動段結(jié)束時刻的高度、速度和彈道傾角.經(jīng)初步仿真發(fā)現(xiàn)，vf主要由助推器的總沖與推力程序決定，受主動段控制規(guī)律影響較?。籋f對總射程影響有限；Θf既對總射程影響大，又對主動段的控制較敏感，可取Θf為主動段彈道與滑翔彈頭飛行彈道之間的銜接量.給定參數(shù)Θf的取值，性能指標(biāo)與彈道約束不變，可得到不同Θf下的彈道曲線以及Θf與全彈道射程Lf的關(guān)系，如圖1、圖2所示.

圖1 不同Θf對應(yīng)的彈道曲線

圖2 全彈道射程與Θf的關(guān)系

可見，Θf對彈道影響較大，且射程隨Θf的變化趨勢并不是簡單的遞增或遞減.滑翔彈頭飛行段彈道呈波狀跳躍形式，當(dāng)Θf小于18°時，全彈道曲線有4個波峰；當(dāng)Θf在18°～25°之間時，彈道曲線有3個波峰；當(dāng)Θf大于25°時，彈道曲線有2個波峰.波峰越多，射程越大.但在以上3種彈道形式中，分別存在3個射程的局部極大值，對應(yīng)的主動段終端傾角為17°，23°，40°，其中Θf＝17°時全彈道射程最大，可將此最優(yōu)傾角記為Θf，opt，其物理意義是保證在同樣的主動段終點機械能條件下，最利于飛行器的能量利用，使射程達到最大值的最優(yōu)傾角.

最大射程優(yōu)化結(jié)果如圖3～圖6所示，可見主動段終點速度略大于4 500m/s，總射程3 300km，滑翔彈頭飛行段射程占總射程的95.5%.落角、落速均嚴格滿足端點約束，法向過載滿足最大過載約束.最大射程彈道呈波狀跳躍形式，大部分時間的飛行高度在30km以上，彈道較平滑.其中主動段結(jié)束至再入大氣層（認為大氣層邊界高度為90km）之前的彈道，認為無控制作用，取控制量為零.

圖3 彈道曲線

圖4 速度曲線

圖5 彈道傾角曲線

圖6 控制量變化曲線

2）起飛質(zhì)量與射程的關(guān)系.

無論是彈道導(dǎo)彈還是助推－滑翔導(dǎo)彈，其全彈道都可分為主動段和被動段，其中被動段又可分為自由飛行段、再入段兩部分.對于彈道導(dǎo)彈而言，再入段射程在整個被動段射程所占比例很小，可近似將該段彈道看成是自由飛行段橢圓彈道的延續(xù)，從而整個被動段射程都用橢圓彈道來計算，其被動段最大射程以及對應(yīng)的最優(yōu)傾角可估算為

式中，vs為第一宇宙速度，R0為平均地球半徑.結(jié)合式（2）、式（15）、式（16）即可得到彈道導(dǎo)彈起飛質(zhì)量與最大射程的關(guān)系.

對于助推－滑翔導(dǎo)彈而言，再入段射程在整個被動段射程中占大部分，而再入段彈道形式復(fù)雜，射程很難用解析式表達.當(dāng)起飛質(zhì)量已知時，取助推火箭燃料比沖Ip＝2 800m/s，火箭結(jié)構(gòu)比ε＝0.1，重推比ν0＝0.5，則可得到助推器基本參數(shù).通過前面的方法，利用彈道優(yōu)化技術(shù)，即得到此起飛質(zhì)量下的最大射程與對應(yīng)的主動段終端傾角.多次計算，即可得出起飛質(zhì)量與最大射程的關(guān)系.計算結(jié)果如圖7～圖10所示.

圖7 不同起飛質(zhì)量下的高度－射程曲線

圖8 不同起飛質(zhì)量下的速度曲線

由圖7、圖8可見，起飛質(zhì)量越大則主動段終點速度越大，全彈道射程也就越大.由圖9可見，在達到同等射程的情況下，助推－滑翔導(dǎo)彈的起飛質(zhì)量要遠比彈道導(dǎo)彈小.由于采用單級火箭助推，在射程大于4 000km以上時，助推－滑翔導(dǎo)彈的起飛質(zhì)量過大，此時宜采用多級火箭助推，以減小消極質(zhì)量.由圖10可見，與彈道導(dǎo)彈相比，助推－滑翔導(dǎo)彈最大射程對應(yīng)的主動段終端傾角要小得多，這就要求導(dǎo)彈主動段具有較強的機動能力，以便完成轉(zhuǎn)彎.

圖9 質(zhì)量比m0/mf隨全彈道射程變化曲線

圖10 最優(yōu)傾角Θf，opt隨m0/mf變化曲線

5 結(jié)論

本文推導(dǎo)了單級固體火箭助推飛行器起飛質(zhì)量與主動段終點速度解析關(guān)系式，給出了助推器參數(shù)的估算方法.使用Gauss偽譜法對助推－滑翔導(dǎo)彈的全彈道進行了優(yōu)化計算.采用分段優(yōu)化的策略，以主動段終端傾角Θf為彈道銜接量，分析了Θf對射程的影響，結(jié)果表明存在使射程達到最大值的最優(yōu)傾角Θf，opt，且此最優(yōu)傾角與彈道導(dǎo)彈相應(yīng)值相比要小得多.最大射程彈道的優(yōu)化結(jié)果表明，Gauss偽譜法處理此類多階段多約束的彈道優(yōu)化問題效果較好，最優(yōu)彈道起伏較小，攻角變化平滑，各項約束都得到滿足.利用彈道優(yōu)化計算，分析了助推－滑翔導(dǎo)彈起飛質(zhì)量與射程的關(guān)系，并與彈道導(dǎo)彈進行了比較，結(jié)果表明助推－滑翔導(dǎo)彈在提高射程方面具有明顯優(yōu)勢.

[1]劉欣，楊濤.滑翔導(dǎo)彈再入拉起段彈道優(yōu)化設(shè)計與制導(dǎo)[J].國防科技大學(xué)學(xué)報，2012，34（1）：67－71.LIU Xin，YANG Tao.Trajectory optimization and guidance in reentry phase for glide missile[J].Journal of National University of Denfence Technology，2012，34（1）：67－71.（in Chinese）

[2]PERSTON H C，DARRY J P.Approximate performance of periodic hypersonic cruise trajectories for global reach[C].Proceedings of AIAA 8th International Space Planes and Hypersonics and Technologies Conference.Norfolk，VA：AIAA，1998：1－6.

[3]李瑜，楊志紅，崔乃剛.助推－滑翔導(dǎo)彈彈道優(yōu)化研究[J].宇航學(xué)報，2008，29（1）：66－71.LI Yu，YANG Zhi－h(huán)ong，CUI Nai－gang.A study of optimal trajectory for boost－glide missle[J].Journal of Astronautics，2008，29（1）：66－71.（in Chinese）

[4]趙吉松，谷良賢，龔春林.高超聲速飛行器起飛質(zhì)量的解析估算[J].固體火箭技術(shù)，2008，31（6）：548－551.ZHAO Ji－song，GU Liang－xian，GONG Chun－lin.Analytical evaluation of takeoff mass for hypersonic vehicles[J].Journal of Solid Rocket Technology，2008，31（6）：548－551.（in Chinese）

[5]賈沛然，陳克俊，何力.遠程火箭彈道學(xué)[M].長沙：國防科技大學(xué)出版社，1993：89－96，195－205.JIA Pei－ran，CHEN Ke－jun，HE Li.Ballistic for long－range rocket[M].Changsha：National University of Denfence Technology Press，1993：89－96，195－205.（in chinese）

[6]TIMOTHY R J.Common aero vehicle autonomous reentry trajectory optimization satisfying waypoint and no－fly zone constraints[D].Ohio，USA：Air Force Insitute of Technology Graduate School of Engineering and Management，2007：88－89.

[7]雍恩米.高超聲速滑翔式再入飛行器軌跡優(yōu)化與制導(dǎo)方法研究[D].長沙：國防科學(xué)技術(shù)大學(xué)，2008：49－53.YONG En－mi.Study on trajectory optimization and guidance approach for hypersonic glide－reentry[D].Changsha：National University of Defense Technology，2008：49－53.（in Chinese）

[8]楊希祥，張為華.基于Gauss偽譜法的固體運載火箭上升段軌跡快速優(yōu)化研究[J].宇航學(xué)報，2011，32（1）：15－21.YANG Xi－xiang，ZHANG Wei－h(huán)ua.Rapid optimization of ascent trajectory for solid launch vehicles based on Gauss pseudospectral method[J].Journal of Astronautics.2011，32（1）：15－21.（in Chinese）

[9]HUNTINGTON G T.Advancement and analysis of a Gauss pseudospectral transcription for optimal control problems[D].Massachusetts，USA：Massachusetts Institute of Technology，2007：51－57，115－143.

[10]BENSON D A，HUNTINGTON G T，THORVALDSEN T P.Direct trajectory optimization and costate estimation via an orthogonal collocation method[J].Journal of Guidance，Control，and Dynamics，2006，29（6）：1 435－1 440.