含滑行時間約束的真空段彈道設(shè)計研究

胡冬生，劉 楠，童科偉，李 爍，張 皓

(1. 中國運載火箭技術(shù)研究院研究發(fā)展中心，北京 100076；2. 中國科學(xué)院空間應(yīng)用工程與技術(shù)中心，北京 100094)

0 引言

在進行運載火箭彈道設(shè)計時，為了提升火箭執(zhí)行太陽同步軌道或近地軌道發(fā)射任務(wù)的性能，通常需要末級具備兩次點火能力，設(shè)置無動力滑行段，以減小速度損失，如我國的長征八號、長征四號丙、長征六號火箭等。未來大量低軌星座任務(wù)的發(fā)射也將依賴于這類彈道設(shè)計。但受到動力系統(tǒng)工程研制、推進劑管理及測控資源等方面因素的制約，末級滑行時間往往受限，無法通過霍曼變軌方式來最大程度發(fā)揮運載火箭的能力[1]。這就帶來了含滑行時間約束的末級多次點火彈道設(shè)計問題，這類問題一般需要火箭滑行時間取約束區(qū)間內(nèi)的最大極限值，以更好地利用重力轉(zhuǎn)彎，因而問題轉(zhuǎn)化為固定滑行時長的彈道設(shè)計。

國內(nèi)在進行帶滑行段飛行任務(wù)的彈道設(shè)計時，通常采用參數(shù)化的工程方法，通過對程序角、關(guān)機時間的迭代和優(yōu)化來獲得滿足入軌約束的最大運載能力[2]。彈道分段增多，設(shè)計變量個數(shù)相應(yīng)增加，使得彈道設(shè)計復(fù)雜，影響了設(shè)計效率。在將最優(yōu)控制理論應(yīng)用于含滑行段的彈道優(yōu)化方面，崔乃剛等[3]研究了大氣層外二三級間有滑行段的最優(yōu)軌跡生成方法，需要對滑行時間進行在線優(yōu)化，無法解決本文提出的問題。

國外Gath等[4]、Lu等[5]、Jezewski[6]主要基于主矢量理論對多次點火飛行器的真空段軌跡和制導(dǎo)設(shè)計進行了研究，但均未對滑行時間進行約束，而是需要優(yōu)化滑行時間，不適用于我國運載火箭設(shè)計的具體情況。Pan等[7]進一步研究推導(dǎo)了含有滑行時間約束的切換條件(Switching Condition)，但未給出應(yīng)用于火箭末級多次點火彈道設(shè)計的方法和說明。

本文在前人研究的基礎(chǔ)上，將含固定滑行時長的切換條件引入真空段含滑行時間約束的末級兩次點火彈道設(shè)計問題中，將彈道優(yōu)化轉(zhuǎn)換成對兩點邊值問題的求解和對運動方程的積分，同時對比分析了兩點邊值問題迭代求解過程中幾種滑行段彈道預(yù)報方法的特點和效果，以及該方法與傳統(tǒng)彈道設(shè)計方法的優(yōu)劣，以期為該類彈道設(shè)計問題提供新的思路和方法。

1 含滑行時間約束的彈道設(shè)計建模

1.1 火箭末級飛行過程

火箭末級一般在大氣層外飛行，通過在兩次工作之間加入滑行段來實現(xiàn)重力轉(zhuǎn)彎，改變火箭飛行方向，從而節(jié)省推進劑。對于兩次點火的火箭來說，要先后經(jīng)歷末級一次動力飛行段、滑行段、末級二次動力飛行段，最終實現(xiàn)入軌，飛行過程見圖1。其中τl1表示末級一次點火時間，由大氣層內(nèi)飛行時間或前一級飛行時間確定；τkl1表示末級一次關(guān)機時間，為設(shè)計變量；τl2表示末級二次點火時間，當(dāng)滑行時間固定時，由τkl1確定；τkl2表示末級二次關(guān)機時間，為設(shè)計變量。

圖1 火箭末級飛行過程Fig.1 Flight process of last stage

1.2 運動方程

火箭在真空段飛行時，所受的力主要是地球引力和推力，其在發(fā)射慣性坐標(biāo)系的運動方程為

(1)

式中，r為火箭相對地心的位置矢量，V為速度矢量；g(r)為引力加速度；T為常值推力；Ib為發(fā)動機推力方向的單位矢量，假定推力方向與火箭縱軸方向重合；m(t)為火箭當(dāng)前質(zhì)量。

采用線性引力場假設(shè)，對運動方程進行歸一化處理[8]，得

(2)

1.3 邊界條件及性能指標(biāo)

(3)

真空段最優(yōu)飛行彈道為推進劑消耗最少的飛行軌跡，在常值推力和滑行時間固定的情況下，即要求飛行時間最短，其性能指標(biāo)設(shè)置如下

minJ=τf

(4)

2 含滑行時間約束的彈道優(yōu)化思路

2.1 最優(yōu)化條件

應(yīng)用最優(yōu)控制理論，哈密頓函數(shù)為

(5)

式中，λr，λV為協(xié)態(tài)變量，ν為拉格朗日乘子。

根據(jù)龐特里亞金極小值原理，最優(yōu)縱軸方向滿足

(6)

協(xié)態(tài)變量滿足

(7)

協(xié)態(tài)變量和狀態(tài)變量的解析解可表示為[9]

(8)

(9)

式中，Ω(τ)，Γ(τ)分別為狀態(tài)轉(zhuǎn)移矩陣、轉(zhuǎn)換矩陣，Ic(τ)，Is(τ)為推力積分，表達式見文獻[9]。當(dāng)火箭處于滑行段時，僅受地球引力作用，Ic(τ)≡0，Is(τ)≡0，可以用二體運動解等彈道預(yù)報方法來代替式(9)，以獲得更為精確的狀態(tài)變量值。

由橫截條件和哈密頓函數(shù)終值條件，最終可得[9]

(10)

(11)

對于滑行時間固定的火箭末級來說，公式(12)的切換條件成立[7]

|λV(τl2)|-|λV(τkl1)|=0

(12)

2.2 程序角計算

在彈道設(shè)計中，一般設(shè)計偏航程序角ψcx=0，軌道傾角約束通過對發(fā)射方位角的迭代來滿足。由此可推導(dǎo)出λr(τ)和λV(τ)的第3個分量均為0，俯仰程序角公式變?yōu)?/p>

(13)

2.3 最優(yōu)彈道求解方法

當(dāng)火箭處于滑行段時，滑行段結(jié)束時刻的位置、速度僅與一次關(guān)機點的位置、速度相關(guān)，因此不進行迭代，僅進行固定滑行時間內(nèi)的彈道積分。

當(dāng)火箭處于末級二次動力飛行段時，通過迭代λr0=[λr0xλr0y0]T，λV0=[λV0xλV0y0]T和τkl25個變量來求解滿足式(3)、(10)、(11)中5個終端約束和最優(yōu)化條件的當(dāng)前協(xié)態(tài)變量、最優(yōu)程序角及關(guān)機時間，再代入原始的運動方程進行彈道積分，進而計算出滿足入軌約束要求和最優(yōu)性的末級二次動力飛行段彈道。該飛行段彈道積分以半長軸約束為終止條件。至此完成了真空段彈道的設(shè)計計算。

具體設(shè)計流程見圖2。圖中，a*表示軌道半長軸約束，a表示當(dāng)前時刻對應(yīng)的半長軸；ε1和ε2為小量，根據(jù)設(shè)計精度要求進行設(shè)定。在各飛行段的迭代計算中，本文采用Levenberg-Marquardt算法進行求解，具有對迭代初值不敏感、求解速度快等優(yōu)點。

圖2 含滑行時間約束的真空段彈道設(shè)計流程Fig.2 Design process of vacuum trajectory with coasting time constraint

3 滑行段彈道預(yù)報方法的影響

由第2章中兩點邊值問題的迭代求解過程可以看出，在末級一次動力飛行段針對每組迭代初值的快速解算過程均需要用到滑行段彈道解析解，本章就幾種預(yù)報方法對滑行段彈道解算和彈道設(shè)計的影響進行分析。除了式(9)所示的解析解外，另外兩種常用的滑行段運動解析解分別為基于二體運動的Goodyear方法[10]和基于中間軌道理論、含J2項攝動的Vinti方法[11-12]。

3.1 滑行段不同彈道預(yù)報方法的精度

若滑行段開始時刻火箭在發(fā)射慣性系的位置為[73 997.5，6 538 780.4，1 303.3]m，速度為[4 080.187，580.715，59.195]m·s-1，則分別利用式(9)解析解、Goodyear方法、含J2項Vinti方法和含J2項數(shù)值積分計算出滑行600 s結(jié)束時的位置、速度，如表1所示。

表1 幾種彈道預(yù)報方法的計算結(jié)果

由表1可見，式(9)解析解與數(shù)值積分的偏差較大，Goodyear結(jié)果偏差較小，Vinti結(jié)果最為接近數(shù)值積分結(jié)果，位置偏差為10 m量級，速度偏差小于0.01 m·s-1。在實際彈道設(shè)計中，可以考慮采用Goodyear或Vinti方法來進行滑行段彈道解算。

3.2 滑行段不同彈道預(yù)報方法對彈道設(shè)計的影響

以文獻[3]中某火箭相關(guān)參數(shù)為算例，輸入?yún)?shù)見表2。末級滑行時間為600 s，目標(biāo)軌道為500 km圓軌道。

表2 真空段彈道設(shè)計輸入?yún)?shù)

分別采用基于Goodyear預(yù)報的彈道設(shè)計方法、基于Vinti預(yù)報的彈道設(shè)計方法開展仿真分析，并與傳統(tǒng)彈道設(shè)計結(jié)果進行對比。計算結(jié)果見表3及圖3～5。

表3 不同彈道設(shè)計方法的計算結(jié)果

圖3 火箭末級飛行地心距Fig.3 Geocentric distance of last stage flight

圖4 火箭末級飛行速度Fig.4 Velocity of last stage flight

圖5 火箭末級飛行俯仰程序角Fig.5 Pitch program angle of last stage flight

在該算例中，傳統(tǒng)彈道設(shè)計方法涉及點火時刻程序角、動力飛行段程序角斜率、關(guān)機時間等6個設(shè)計變量，約束參數(shù)僅有入軌點地心距、速度和彈道傾角3個，遠少于設(shè)計變量個數(shù)，因此需要進行大量、充分的優(yōu)化運算；而本文提出的方法僅需針對兩點邊值問題進行不斷的迭代計算和彈道積分，思路相對簡潔。仿真結(jié)果表明，3種方法均可以得到滿足入軌約束的飛行彈道，關(guān)機時間、入軌質(zhì)量較為接近，其中基于Vinti的彈道設(shè)計結(jié)果與傳統(tǒng)方法設(shè)計結(jié)果更為吻合，包括一次、二次關(guān)機時間等，最終的入軌質(zhì)量僅相差8.7 kg。而基于Goodyear的彈道設(shè)計方法采用了線性引力場假設(shè)和滑行段二體運動解析解，迭代模型精度有所損失，進而導(dǎo)致了在偏心率的入軌約束上精度稍差，入軌質(zhì)量較其他2種方法小大約180 kg，二次動力飛行段程序角斜率也更大，達到約8(°)/s， 將給火箭姿態(tài)控制帶來不利影響。

4 多級火箭含滑行時間約束的彈道設(shè)計

4.1 多級火箭含滑行時間約束的彈道設(shè)計思路

對于末級可以兩次啟動的三級火箭而言，通常二級和三級均在大氣層外飛行，本文提出的彈道設(shè)計方法仍然適用，只需在2.3節(jié)的設(shè)計流程中增加二級飛行段。由于二級飛行時間固定，因此迭代參數(shù)不變，通過設(shè)置迭代初值λr0=[λr0xλr0y0]T，λV0=[λV0xλV0y0]T,τkl1，τkl2和式(8)、(9)依次快速解算出二級飛行段、三級一次動力飛行段、滑行段、三級二次動力飛行段的狀態(tài)變量和協(xié)態(tài)變量，其余過程與2.3節(jié)相同。

本節(jié)以某常溫三級火箭為例進行仿真計算[13]。

4.2 彈道設(shè)計算例仿真

火箭二級和三級相關(guān)參數(shù)取值見表4。三級滑行時間為550 s，目標(biāo)軌道為600 km太陽同步軌道。

分別采用基于Goodyear預(yù)報的彈道設(shè)計方法、基于Vinti預(yù)報的彈道設(shè)計方法、傳統(tǒng)彈道設(shè)計方法開展仿真分析。計算結(jié)果見表5及圖6～8。

表4 三級火箭真空段彈道設(shè)計輸入?yún)?shù)

表5 三級火箭不同彈道設(shè)計方法的計算結(jié)果

圖6 火箭二、三級飛行地心距Fig.6 Geocentric distance of second and third stage flight

圖7 火箭二、三級飛行速度Fig.7 Velocity of second and third stage flight

圖8 火箭二、三級飛行俯仰程序角Fig.8 Pitch program angle of second and third stage flight

在該算例中，采用傳統(tǒng)彈道設(shè)計方法時由于增加了二級飛行段，設(shè)計變量增加到8個，而約束參數(shù)仍然只有入軌點地心距、速度和彈道傾角3個，這就進一步加大了優(yōu)化難度；而本文提出的方法相比3.2節(jié)算例并沒有增加迭代變量的個數(shù)，求解過程幾乎完全一致，僅增加了二級飛行段的迭代和彈道積分。

仿真結(jié)果表明，基于Vinti的彈道設(shè)計結(jié)果與傳統(tǒng)方法設(shè)計結(jié)果更為吻合，最終的入軌質(zhì)量僅相差0.8 kg；而基于Goodyear的彈道設(shè)計結(jié)果在入軌質(zhì)量上較其他兩種方法減小近70 kg，二次動力飛行段程序角斜率也更大，達到3.09(°)/s，與3.2節(jié)算例出現(xiàn)的情況一致。

5 兩種彈道設(shè)計思路對比

傳統(tǒng)彈道設(shè)計方法本質(zhì)上是將程序角、關(guān)機時間等設(shè)計變量代入完整、準(zhǔn)確的運動方程進行數(shù)值積分，進而開展迭代和優(yōu)化，該方法中的程序角采用了工程上常用的線性化處理，在程序角設(shè)計規(guī)律上有時會損失一定的最優(yōu)性。此外，彈道設(shè)計的最優(yōu)性主要靠優(yōu)化算法來保證，優(yōu)化算法的魯棒性、尋優(yōu)能力都將影響最終的彈道設(shè)計結(jié)果。尤其是在含無動力滑行段的情況下，設(shè)計變量個數(shù)遠遠多于約束個數(shù)，滑行段的存在也進一步增強了問題的非線性，使得彈道優(yōu)化難度加大，對程序角初值的設(shè)置更為敏感，甚至陷入局部最優(yōu)解。

本文提出的基于線性引力場和滑行段切換條件的彈道設(shè)計方法相對簡潔很多，其基本的迭代變量僅有初始的4個協(xié)態(tài)初值和2個關(guān)機時間，彈道設(shè)計的最優(yōu)性完全依靠橫截條件、哈密頓函數(shù)終值條件和切換條件來保證。但同時，該方法在迭代過程中使用了簡化的運動方程解析解，且最優(yōu)性條件也是基于線性引力場推導(dǎo)出來的，這也在一定程度上損失了最優(yōu)性。此外，通過仿真可以發(fā)現(xiàn)，滑行段彈道預(yù)報模型越精確，則迭代產(chǎn)生的程序角越接近最優(yōu)解，最終得到的入軌質(zhì)量就越大，也越容易滿足對入軌軌道的約束精度，同時也不會產(chǎn)生較大的程序角斜率。因此，在含滑行段的彈道設(shè)計中，應(yīng)當(dāng)優(yōu)先采用Vinti預(yù)報方法開展滑行段彈道解算。

通過對比和仿真可知，兩種彈道設(shè)計方法各有特點，通常可以獲得幾乎相同的設(shè)計結(jié)果。當(dāng)應(yīng)用于級數(shù)較多、含滑行段的火箭彈道設(shè)計問題時，基于線性引力場和切換條件的彈道設(shè)計方法能夠在一定程度上簡化問題的求解，保證得到問題的近似最優(yōu)解，且在小推重比情況下更具有性能優(yōu)勢[8]；此外，使用的Vinti滑行段彈道預(yù)報方法也可以應(yīng)用在傳統(tǒng)彈道設(shè)計的滑行段計算中，從而避免滑行段彈道積分，大大提高彈道設(shè)計的效率。

6 結(jié)論

本文在線性引力場假設(shè)的基礎(chǔ)上，將含固定滑行時長的切換條件引入彈道設(shè)計中，從而將含滑行時間約束的彈道優(yōu)化問題轉(zhuǎn)換成對兩點邊值問題的求解和對運動方程的積分。通過多個算例仿真得出以下結(jié)論：

1)基于線性引力場和滑行段切換條件的彈道設(shè)計方法思路簡潔，能夠適用于對滑行時間有約束、末級兩次點火的兩級或多級火箭，設(shè)計結(jié)果與傳統(tǒng)設(shè)計方法幾乎一致；

2)滑行段彈道預(yù)報精度對最終的彈道優(yōu)化結(jié)果具有一定的影響，模型精度越高，優(yōu)化結(jié)果越好；

3)本文方法和思路應(yīng)用范圍較廣，能夠在各種推重比情況下獲得最優(yōu)彈道，在總體方案論證或初步設(shè)計階段可以替代傳統(tǒng)設(shè)計方法，以保證獲得更好的運載性能。