一種針對(duì)直接入軌飛行器的主動(dòng)段軌道快速實(shí)現(xiàn)方案

王 鵬,李 強(qiáng),王 智,宋劍爽,寧 雷

(北京宇航系統(tǒng)工程研究所，北京 100076)

0 引言

軌道設(shè)計(jì)作為飛行器發(fā)射中一項(xiàng)重要的技術(shù)保障，直接關(guān)系到飛行器的發(fā)射準(zhǔn)備時(shí)間和性能指標(biāo)，進(jìn)而影響到執(zhí)行任務(wù)的能力?？紤]性能指標(biāo)的主動(dòng)段軌道設(shè)計(jì)屬于帶過(guò)程約束的終端時(shí)刻自由、終端狀態(tài)受約束的非線性最優(yōu)控制問(wèn)題，早期多采用間接法將最優(yōu)控制問(wèn)題轉(zhuǎn)化為邊值問(wèn)題，進(jìn)而通過(guò)數(shù)值方法求解，隨著計(jì)算機(jī)水平的提高，目前多采用直接法求解，在不考慮時(shí)間指標(biāo)的情況下均取得了較好的使用效果[1]。

為提升飛行器射前快速反應(yīng)能力，希望在保證一定性能指標(biāo)的情況下盡可能縮短軌道計(jì)算時(shí)間。諸多學(xué)者針對(duì)此問(wèn)題進(jìn)行了相關(guān)研究，當(dāng)前主流的方法是通過(guò)離線計(jì)算或全局優(yōu)化算法獲得合理的參數(shù)初值，以此為基礎(chǔ)采用序列二次規(guī)劃等梯度數(shù)值算法求解[2]。考慮到相關(guān)方法對(duì)于優(yōu)化指標(biāo)并沒(méi)有放松，因此仍屬于最優(yōu)化問(wèn)題，所耗費(fèi)的時(shí)間仍然比較長(zhǎng)。

本文以直接入軌飛行器為研究對(duì)象，在不考慮地球旋轉(zhuǎn)的情況下通過(guò)優(yōu)化算法求解得到具有優(yōu)化指標(biāo)的軌道設(shè)計(jì)結(jié)果，在旋轉(zhuǎn)橢球模型下基于球面三角確定射向，以此為基礎(chǔ)考慮旋轉(zhuǎn)影響對(duì)地球不轉(zhuǎn)情況下的設(shè)計(jì)結(jié)果進(jìn)行修正得到次優(yōu)化的軌道結(jié)果。由于優(yōu)化過(guò)程可以在射前離線進(jìn)行，且射前僅通過(guò)參數(shù)修正即可得到當(dāng)前射向?qū)?yīng)的軌道設(shè)計(jì)參數(shù)，省去了射前的軌道迭代計(jì)算過(guò)程，因此發(fā)射準(zhǔn)備時(shí)間大大縮短，有效地提升了射前快速反應(yīng)能力。

1 問(wèn)題描述

1.1 動(dòng)力學(xué)方程

忽略內(nèi)部介質(zhì)相對(duì)于飛行器流動(dòng)所引起的附加科氏力，假設(shè)在控制系統(tǒng)作用下飛行器始終處于力矩瞬時(shí)平衡狀態(tài)，忽略控制力的影響，慣性坐標(biāo)系中以矢量描述的質(zhì)心動(dòng)力學(xué)方程為

(1)

式中，P為推力矢量，R為氣動(dòng)力矢量，mg為引力矢量。設(shè)動(dòng)參考系相對(duì)于慣性坐標(biāo)系以角速度ωe轉(zhuǎn)動(dòng)，由矢量導(dǎo)數(shù)法可知

(2)

將式(1)代入式(2)，可得動(dòng)參考系中的質(zhì)心運(yùn)動(dòng)學(xué)方程[3]

(3)

將上述各項(xiàng)投影到發(fā)射坐標(biāo)系中即可得到發(fā)射坐標(biāo)系中的動(dòng)力學(xué)方程。

1.2 設(shè)計(jì)變量選取

直接入軌飛行器主動(dòng)段軌道設(shè)計(jì)與傳統(tǒng)軌道飛行器相似，可采用相同的飛行程序角模型，如圖1所示。

圖1 主動(dòng)段飛行程序角Fig.1 Boost phase flight program angle

即垂直起飛后不久開始在稠密大氣層內(nèi)以αmax實(shí)行負(fù)攻角轉(zhuǎn)彎，在飛行速度接近跨聲速段之前，將攻角收縮為零，在重力作用下繼續(xù)轉(zhuǎn)彎，直到滿足分離條件φ3，為提高飛行器的穩(wěn)定性，在發(fā)動(dòng)機(jī)啟動(dòng)、關(guān)機(jī)以及級(jí)間分離前后，要求飛行器盡量減小飛行攻角，分離后通過(guò)末級(jí)飛行段設(shè)計(jì)攻角αm1和αm2進(jìn)行能量管理以滿足交班點(diǎn)約束條件。據(jù)此，飛行程序角模型可取為

φ=φ(αmax,φ3,αm1,αm2)

(4)

式中，φ為俯仰程序角，因[t2,t3]段重力轉(zhuǎn)彎無(wú)設(shè)計(jì)變量，而在無(wú)動(dòng)力滑行時(shí)，φ3和t3取一個(gè)即可，這里將φ3作為設(shè)計(jì)變量。同時(shí)，由于地球形狀不規(guī)則以及自轉(zhuǎn)影響，主動(dòng)段軌道還與發(fā)射點(diǎn)地理緯度B0、高度h0和發(fā)射方位角A0有關(guān)，其中發(fā)射點(diǎn)是已知的，發(fā)射方位角可根據(jù)規(guī)劃要求計(jì)算得到。

1.3 約束條件

主動(dòng)段約束條件包括：

1)飛行最大動(dòng)壓約束：

Q≤Qmax

(5)

2)級(jí)間分離起控條件：

(6)

3)終端交班點(diǎn)高度和當(dāng)?shù)剀壍纼A角約束：

|hf-hfstd|≤Δhf，|Θf|≤ΔΘf

(7)

其中，Qmax為飛行最大動(dòng)壓約束；下標(biāo)fl表示分離點(diǎn)參數(shù)，hflmin為級(jí)間分離最低高度，Qflmax為級(jí)間分離最大動(dòng)壓；下標(biāo)f表示終端參數(shù)，Δhf為終端高度偏差，ΔΘf為終端當(dāng)?shù)剀壍纼A角約束。

1.4 優(yōu)化指標(biāo)

優(yōu)化指標(biāo)為終端的速度最大，即

J=-Vf

(8)

2 軌道優(yōu)化設(shè)計(jì)與實(shí)現(xiàn)

在主動(dòng)段滿足飛行約束助推至交班點(diǎn)后，考慮到直接入軌飛行器自身可在入軌后大范圍機(jī)動(dòng)飛行，因此交班點(diǎn)的條件并非像傳統(tǒng)衛(wèi)星軌道要求那樣嚴(yán)苛，在交班點(diǎn)高度和當(dāng)?shù)貜椀纼A角適當(dāng)放寬的情況下，可轉(zhuǎn)化為滿足約束要求的可行區(qū)域。交班點(diǎn)條件的放寬帶來(lái)了主動(dòng)段軌道實(shí)現(xiàn)的新思路，即在不考慮地球自轉(zhuǎn)、軌道設(shè)計(jì)與射向無(wú)關(guān)的情況下進(jìn)行軌道相關(guān)參數(shù)的優(yōu)化設(shè)計(jì)，實(shí)際使用過(guò)程中基于射向和發(fā)射點(diǎn)參數(shù)對(duì)軌道設(shè)計(jì)參數(shù)進(jìn)行一定的修正，以降低交班點(diǎn)參數(shù)的偏差，最終滿足交班點(diǎn)的約束條件。

2.1 不考慮地球旋轉(zhuǎn)的設(shè)計(jì)參數(shù)優(yōu)化

本文研究的內(nèi)容為若干個(gè)設(shè)計(jì)變量情況下滿足約束條件的參數(shù)優(yōu)化問(wèn)題，優(yōu)化算法可選用序列二次規(guī)劃算法(SQP)，在離線優(yōu)化過(guò)程中參數(shù)初值可通過(guò)經(jīng)驗(yàn)法選取。

SQP算法的基本思想是在給定的近似點(diǎn)處通過(guò)二次近似逐漸得到一個(gè)更好的迭代點(diǎn)，這需要通過(guò)求解二次規(guī)劃子問(wèn)題得到，在當(dāng)前迭代點(diǎn)處算法通過(guò)求解一系列的二次規(guī)劃子問(wèn)題，使得迭代點(diǎn)逐步接近原優(yōu)化命題的最優(yōu)點(diǎn)，算法最終收斂到最優(yōu)解[4]。

設(shè)f,g,h分別為目標(biāo)函數(shù)、等式約束和不等式約束且均二次連續(xù)可微，非線性規(guī)劃問(wèn)題(P)可表示為

(9)

對(duì)于問(wèn)題(P)，SQP方法通過(guò)序列地求解一系列的二次規(guī)劃子問(wèn)題來(lái)逐步逼近原問(wèn)題的最優(yōu)解，在迭代點(diǎn)xk處，與上式相對(duì)應(yīng)的二次規(guī)劃子問(wèn)題(QP)可表示為

(10)

為實(shí)現(xiàn)SQP算法，還須解決兩個(gè)問(wèn)題[5]，其一是矩陣Bk的修正，其二是QP子問(wèn)題相容性。

在不考慮地球旋轉(zhuǎn)的情況下，軌道優(yōu)化設(shè)計(jì)過(guò)程中無(wú)須考慮發(fā)射點(diǎn)經(jīng)緯度和射向的影響，由此可優(yōu)化得到給定發(fā)射點(diǎn)高度下的軌道設(shè)計(jì)結(jié)果，作為考慮旋轉(zhuǎn)橢球模型下設(shè)計(jì)參數(shù)修正的輸入。

2.2 考慮地球旋轉(zhuǎn)的設(shè)計(jì)參數(shù)修正

圖2 飛行過(guò)程中的受力情況Fig.2 The force of flight process

對(duì)于直接入軌飛行器而言，相關(guān)的控制指令主要是俯仰程序角，在不考慮三通道參數(shù)相互影響的情況下，對(duì)俯仰程序角和飛行過(guò)程中的力對(duì)軌道的影響進(jìn)行定性分析：

1)在同一套發(fā)射慣性系俯仰角φT(t)控制參數(shù)下，進(jìn)行發(fā)射系俯仰角變化分析：

① 當(dāng)A0∈(0°,180°)，即射向向東時(shí)，任意時(shí)刻的φ(t)>φT(t)；

② 當(dāng)A0∈(180°,360°)，即射向向西時(shí)，任意時(shí)刻的φ(t)<φT(t)；

③ 當(dāng)A0=0°或180°，即射向北或向南時(shí)，任意時(shí)刻的φ(t)=φT(t)。

2)在同一套發(fā)射系俯仰角φ(t)控制參數(shù)下，進(jìn)行受力分析?？紤]到在接近交班點(diǎn)過(guò)程中，速度方向逐漸接近當(dāng)?shù)厮骄€，此時(shí)：

① 當(dāng)V向東時(shí)，科氏慣性力方向向上，等效于減小引力加速度的作用，交班點(diǎn)高度偏高；

② 當(dāng)V向西時(shí)，科氏慣性力方向向下，等效于增加引力加速度的作用，交班點(diǎn)高度偏低；

③ 當(dāng)V向北或南時(shí)，科氏慣性力為0，交班點(diǎn)高度保持不變。

由上述分析可知，對(duì)于同一套發(fā)慣系俯仰角φT(t)，受地球旋轉(zhuǎn)的影響，不但實(shí)際的控制指令由φT(t)變?yōu)棣?t)，而且即使對(duì)于同一套φ(t)，受附加力的影響，軌道的終端高度也會(huì)發(fā)生變化。

上述兩方面影響中，程序角可以進(jìn)行人為修正，由φ(t)=φT(t)-ωe(-cosB0sinA0)t可知，令

φ′T(t)=φT(t)+ωe(-cosB0sinA0)t

(11)

即可消除程序角指令變化帶來(lái)的影響，其中φ′T(t)為修正后的發(fā)射慣性系俯仰角。

受力關(guān)系雖然不可補(bǔ)償，但注意到上述程序角指令變化和產(chǎn)生的附加力對(duì)終端高度的影響都是同向的，例如當(dāng)向東發(fā)射時(shí)，任意時(shí)刻φ(t)>φT(t)本身已經(jīng)使得終端高度變高，而科氏慣性力的影響則使得終端的高度變得更高，因此可通過(guò)進(jìn)一步修正程序角參數(shù)來(lái)“補(bǔ)償”地球自轉(zhuǎn)引起的附加力作用下的軌道偏差。

綜上所述，修正項(xiàng)可選為K·ωe(-cosB0sinA0)t，則修正后的發(fā)射慣性系系程序角為

φ′T(t)=φT(t)+Kωe(-cosB0sinA0)t(K>1)

(12)

考慮到上述修正項(xiàng)中已經(jīng)包含了時(shí)間t，可消除時(shí)間增長(zhǎng)對(duì)后續(xù)軌道的影響，因此K一般取常值即可。

2.3 主動(dòng)段軌道實(shí)現(xiàn)方案

考慮發(fā)射點(diǎn)高度的影響，通過(guò)在不考慮地球旋轉(zhuǎn)的情況下進(jìn)行軌道優(yōu)化設(shè)計(jì)得到相關(guān)設(shè)計(jì)結(jié)果，將程序角參數(shù)考慮地球旋轉(zhuǎn)影響進(jìn)行修正，則可得到適用于旋轉(zhuǎn)橢球下的程序角參數(shù)，由此可形成主動(dòng)段軌道實(shí)現(xiàn)方案：

1)在地球不轉(zhuǎn)情況下，考慮發(fā)射點(diǎn)高度區(qū)間，按高度間隔分別使用SQP優(yōu)化算法計(jì)算得到發(fā)射點(diǎn)高度對(duì)應(yīng)的程序角變化參數(shù)，最終形成以發(fā)射點(diǎn)高度為自變量的一維程序角數(shù)組；

2)基于實(shí)際的發(fā)射點(diǎn)高度插值程序角數(shù)組得到此高度對(duì)應(yīng)的地球不轉(zhuǎn)情況下的程序角參數(shù)；

3)考慮直接入軌飛行器機(jī)動(dòng)飛行要求，根據(jù)交班點(diǎn)經(jīng)緯度信息，采用球面三角理論確定發(fā)射方位角；

4)考慮地球旋轉(zhuǎn)影響，基于發(fā)射方位角和發(fā)射點(diǎn)緯度，對(duì)程序角參數(shù)進(jìn)行修正，得到最終程序角參數(shù)。整個(gè)主動(dòng)段軌道實(shí)現(xiàn)方案的執(zhí)行流程如圖3所示。

圖3 主動(dòng)段軌道實(shí)現(xiàn)方案Fig.3 Trajectory implementation of boost phase

3 仿真算例與分析

為驗(yàn)證上述軌道實(shí)現(xiàn)方案的有效性和可行性，基于地球不轉(zhuǎn)模型，在發(fā)射點(diǎn)經(jīng)緯度均為0，初始高度為h0、發(fā)射點(diǎn)高度間隔取h1的情況下，利用優(yōu)化算法分別得到滿足交班點(diǎn)參數(shù)要求的不同發(fā)射點(diǎn)高度情況下隨相對(duì)發(fā)射點(diǎn)高度變化的發(fā)慣系俯仰角程序角插值數(shù)組，見(jiàn)圖4。交班點(diǎn)參數(shù)要求包括終端高度和終端當(dāng)?shù)貜椀纼A角兩個(gè)要求：其中終端當(dāng)?shù)貜椀纼A角要求為0°；關(guān)于終端高度，優(yōu)化過(guò)程中要求嚴(yán)格等于標(biāo)準(zhǔn)值，而在程序角修正后的終端高度容許在標(biāo)準(zhǔn)值的基礎(chǔ)上存在3 km的高度偏差。在上述基礎(chǔ)上，交班點(diǎn)的速度越大越好。

基于上述程序角插值數(shù)組，已知發(fā)射點(diǎn)高度參數(shù)，可插值獲取當(dāng)前高度下不考慮地球旋轉(zhuǎn)的程序角參數(shù)，采用上述程序角修正方案對(duì)相關(guān)程序角進(jìn)行修正即可得到當(dāng)前高度下考慮地球旋轉(zhuǎn)的實(shí)際程序角參數(shù)。

圖4 發(fā)慣系俯仰角隨相對(duì)發(fā)射點(diǎn)高度變化曲線Fig.4 Relationship between flight program angle and height

3.1 計(jì)算時(shí)間對(duì)比

在Windows XP操作系統(tǒng)上，采用主頻≥2.6 GHz的I5處理器、容量≥4 GB的DDR2/3內(nèi)存，基于C語(yǔ)言編寫的計(jì)算程序，在同一發(fā)射點(diǎn)和交班狀態(tài)參數(shù)下，采用程序角迭代計(jì)算和本文所使用的程序角插值后修正分別計(jì)算一條軌道所需的最長(zhǎng)時(shí)間,對(duì)比結(jié)果如表1所示。

表1 計(jì)算時(shí)間對(duì)比

根據(jù)上述對(duì)比可知，采用程序角插值修正的方式由于沒(méi)有迭代計(jì)算過(guò)程，大大縮短了計(jì)算量，使得計(jì)算一條軌道所需的時(shí)間顯著減少，能夠?qū)崿F(xiàn)快速計(jì)算的要求。

3.2 交班點(diǎn)偏差分布對(duì)比

采用模擬打靶方式對(duì)國(guó)內(nèi)緯度[18°，53°]和高度[0 m,3 000 m]范圍內(nèi)任意發(fā)射點(diǎn)、任意射向情況下進(jìn)行軌道計(jì)算，計(jì)算過(guò)程中未考慮制導(dǎo)修正，以最高點(diǎn)作為終止計(jì)算條件，直接使用地球不轉(zhuǎn)情況下的程序角和修正后的程序角分別計(jì)算得到的交班點(diǎn)高度偏差和速度偏差分布見(jiàn)圖5和圖6，相對(duì)標(biāo)準(zhǔn)值的標(biāo)準(zhǔn)差計(jì)算結(jié)果見(jiàn)表2，其中標(biāo)準(zhǔn)值是指不考慮地球旋轉(zhuǎn)情況下通過(guò)優(yōu)化得到的滿足交班點(diǎn)約束的數(shù)值。

圖5 交班點(diǎn)高度偏差散布圖Fig.5 Intersection height scatter

圖6 交班點(diǎn)速度偏差散布圖Fig.6 Intersection speed scatter

表2 標(biāo)準(zhǔn)差結(jié)果

由以上結(jié)果可得，相對(duì)不考慮地球旋轉(zhuǎn)情況下的標(biāo)準(zhǔn)值，受點(diǎn)位和射向影響考慮地球旋轉(zhuǎn)情況下速度和高度都會(huì)在標(biāo)準(zhǔn)值附近出現(xiàn)散布，但使用修正后的程序角參數(shù)可顯著降低地球旋轉(zhuǎn)因素的影響，最終交班點(diǎn)的高度和速度偏差大幅降低，在滿足交班點(diǎn)要求的同時(shí)為后續(xù)飛行創(chuàng)造良好的初始條件。

3.3 優(yōu)化指標(biāo)情況對(duì)比

為驗(yàn)證修正程序角后軌道計(jì)算結(jié)果的優(yōu)化指標(biāo)情況，選取h0高度下的程序角參數(shù)，采用前述方法修正程序角，計(jì)算發(fā)射方位角A0∈[0,360)范圍內(nèi)的交班點(diǎn)參數(shù)。同時(shí)以此發(fā)射點(diǎn)參數(shù)和發(fā)射方位角作為初始條件，考慮過(guò)程約束和終端交班點(diǎn)約束，基于優(yōu)化算法計(jì)算同樣發(fā)射方位角范圍下的交班點(diǎn)參數(shù)，兩者高度偏差和速度偏差的對(duì)比結(jié)果見(jiàn)圖7和圖8。

圖7 交班點(diǎn)高度偏差隨發(fā)射方位角變化情況Fig.7 Relationship between intersection height and launching direction

圖8 交班點(diǎn)速度偏差隨發(fā)射方位角變化情況Fig.8 Relationship between intersection speed and launching direction

由圖7、圖8可以看到，采用SQP算法優(yōu)化可以精確滿足終端高度約束條件，因此不隨發(fā)射方位角變化，而終端速度受地球自轉(zhuǎn)影響隨發(fā)射方位角先增大后減小。對(duì)于采用程序修正的結(jié)果，終端高度和速度受能量守恒約束呈現(xiàn)此消彼長(zhǎng)的規(guī)律，在不考慮終端高度差異的情況下，較優(yōu)化結(jié)果使用程序角修正的最大速度偏小2.5 m/s左右，若考慮高度差異，則兩者的速度差異將更小。由此說(shuō)明，采用程序角修正的軌道實(shí)現(xiàn)方法仍具有一定的優(yōu)化指標(biāo)。

4 結(jié)論

本文研究了一種針對(duì)直接入軌飛行器的主動(dòng)段軌道實(shí)現(xiàn)方案，采用地球不轉(zhuǎn)模型計(jì)算指定高度下的程序角參數(shù)，使用時(shí)根據(jù)發(fā)射點(diǎn)高度插值得到對(duì)應(yīng)的程序角，基于發(fā)射點(diǎn)緯度和射向?qū)Υ顺绦蚪沁M(jìn)行修正后作為考慮地球旋轉(zhuǎn)模型的實(shí)際程序角。仿真結(jié)果表明，該方法僅需進(jìn)行離線軌道優(yōu)化，無(wú)需進(jìn)行實(shí)時(shí)迭代計(jì)算，大大縮短了軌道實(shí)現(xiàn)所需的準(zhǔn)備時(shí)間，同時(shí)所計(jì)算出的交班點(diǎn)參數(shù)偏差較小，能夠滿足交班點(diǎn)要求且具有一定的優(yōu)化指標(biāo)，具有一定的工程參考價(jià)值。