王 鵬,李 強,王 智,宋劍爽,寧 雷

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

0 引言

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

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

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

1 問題描述

1.1 動力學(xué)方程

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

(1)

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

(2)

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

(3)

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

1.2 設(shè)計變量選取

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

圖1 主動段飛行程序角Fig.1 Boost phase flight program angle

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

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

(4)

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

1.3 約束條件

主動段約束條件包括：

1)飛行最大動壓約束：

Q≤Qmax

(5)

2)級間分離起控條件：

(6)

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

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

(7)

其中，Qmax為飛行最大動壓約束；下標(biāo)fl表示分離點參數(shù)，hflmin為級間分離最低高度，Qflmax為級間分離最大動壓；下標(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è)計與實現(xiàn)

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

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

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

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

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

(9)

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

(10)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(11)

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

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

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

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

(12)

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

2.3 主動段軌道實現(xiàn)方案

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

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

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

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

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

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

3 仿真算例與分析

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

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

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

3.1 計算時間對比

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

表1 計算時間對比

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

3.2 交班點偏差分布對比

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

圖5 交班點高度偏差散布圖Fig.5 Intersection height scatter

圖6 交班點速度偏差散布圖Fig.6 Intersection speed scatter

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

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

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

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

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

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

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

4 結(jié)論

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