助推－滑翔導(dǎo)彈彈道優(yōu)化與總體參數(shù)分析

劉 欣，楊 濤，張青斌

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3）約束條件.

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

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

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

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

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

1）時(shí)域變換.

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

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

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

Gauss偽譜法的離散點(diǎn)κ＝｛τ1，…，τK｝為K階Legendre－Gauss點(diǎn)，即K階Legendre多項(xiàng)式的根PK（τ），其中，

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

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

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

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

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

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

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

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

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

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

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

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

則連接點(diǎn)約束條件為

基于上述的數(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)性條件之間的等價(jià)性.除此之外，GPM的主要優(yōu)點(diǎn)是算法以指數(shù)速率收斂，收斂速率快，且能在離散點(diǎn)較少時(shí)得到較好的精度.

4 優(yōu)化算例

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

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

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

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

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

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

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

圖3 彈道曲線

圖4 速度曲線

圖5 彈道傾角曲線

圖6 控制量變化曲線

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

[1]劉欣，楊濤.滑翔導(dǎo)彈再入拉起段彈道優(yōu)化設(shè)計(jì)與制導(dǎo)[J].國(guó)防科技大學(xué)學(xué)報(bào)，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é)報(bào)，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]賈沛然，陳克俊，何力.遠(yuǎn)程火箭彈道學(xué)[M].長(zhǎng)沙：國(guó)防科技大學(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].長(zhǎng)沙：國(guó)防科學(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ùn)載火箭上升段軌跡快速優(yōu)化研究[J].宇航學(xué)報(bào)，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.