一種常值推力作用的軌道機動優化設計策略

劉將輝，彭祺擘，張海聯，周建平

（1.中國航天員科研訓練中心，北京 100094；2.中國載人航天工程辦公室，北京 100071；3.北京航空航天大學宇航學院，北京 100083）

0 引言

近年來，隨著航天應用技術和深空探測的迅速發展，作為實現空間復雜任務的基礎，連續推力作用下的航天器軌道機動優化設計問題在空間探測中具有越來越重要的意義［1-3］。

針對連續推力的航天器軌道機動問題，相關研究主要集中在常值推力和非常值推力這兩種不同類型推力作用下的軌道機動方式。常值推力作用的軌道機動指航天器在大小恒定、方向可控的推力作用下實現軌道轉移的機動方式［4-6］。文獻［7］研究了連續常值推力在軌道機動中的應用，該研究的軌道機動方式僅限于共面軌道轉移和共面軌道逃逸。文獻［8］研究了二體條件下連續常值徑向推力作用的軌道機動問題，得到了航天器在常值徑向推力作用時沿圓軌道飛行的條件。文獻［9］研究了多段常值推力的水滴懸停軌道控制問題。文獻［10］研究了常值推力作用的航天器最優交會問題，通過分段多項式描述航天器的狀態量和控制量。文獻［11］采用間接法研究了連續常值推力作用的共面圓軌道最小時間軌道轉移問題。文獻［12］采用改進間接法研究了常值推力作用下航天器共面軌道轉移的燃料優化問題，給出了最優軌跡的方向角約束條件。文獻［13］分析了航天器在常值徑向推力作用下的軌道運動特性。文獻［14］研究了考慮通信窗口約束的常值推力橢圓軌道最優交會問題。目前有關常值推力作用的軌道研究雖然取得許多成果，但局限于面內軌道機動，使得其存在應用范圍過窄的缺陷。

非常值推力作用的軌道機動研究方面，文獻［15］研究了連續推力作用的人工太陽同步軌道設計問題。文獻［16］分析了連續小推力作用下航天器軌道機動動力學特性。文獻［17-18］將軌道轉移問題轉化成凸優化問題，提出了一種分布式空間系統自主重構的新方法，該方法可確保航天器編隊的安全性，同時可優化燃料消耗。針對近地圓軌道航天器編隊重構的軌道機動問題，文獻［19］提出了考慮速度增量約束的多脈沖機動規劃算法，文獻［20］提出一種基于凸優化方法的最優軌跡規劃方法。文獻［21-22］研究了連續小推力軌道機動方式在空間碎片清除任務中的應用。文獻［23］分析了將大型空間碎片移軌/轉軌到處置軌道的兩種策略效率。文獻［24］研究了面向小行星探測的連續小推力變軌問題，提出了采用自適應交互式多模型無跡卡爾曼濾波算法。文獻［25］研究了微推力對微納衛星軌道機動的影響，文獻［26］分析了小推力和大氣阻力對近地軌道機動的影響。文獻［27］研究了圓軌道小推力衛星的軌道機動與燃料優化問題，文獻［28］研究了采用電推進的小推力地心軌道轉移的燃料優化問題。針對靜止軌道衛星的軌道機動問題，文獻［29］提出一種4 脈沖遠程快速交會規劃方法。文獻［30］研究了航天器在近圓軌道之間機動的混合推力問題，提出了一種改進的交互多模型算法。文獻［31］首次提出基于逆向設計的形狀法，假設航天器沿某一類型的曲線飛行，通過曲線擬合法求解方程中的待定系數，從而得到航天器飛行軌跡的解析表達式，實現對連續小推力航天器軌道機動問題的快速計算。文獻［32］采用正弦指數模型設計了小推力星際轉移軌道。針對帶推力約束的航天器空間軌道機動問題，文獻［33］提出了有限傅里葉級數模型，文獻［34］提出了一種改進的形狀設計方法。文獻［35］將軌道形狀方程與貝塞爾曲線方程應用到平面軌道機動問題中，提出了一種基于貝塞爾曲線的形狀設計方法。

相比于常值推力作用的軌道研究，非常值推力研究拓展到三維空間的軌道機動，但存在3點不足：1）工程實際中航天器的推力大小一般是恒定的。2）空間軌道機動的研究大多采用脈沖變軌這種理想化的變軌方式，該方式僅適用于機動時長占總飛行時長比例較小的情況。3）目前非常值推力的研究大多采用基于逆向設計的形狀法，該方法適合于面內圓軌道之間的軌道轉移，當軌道偏心率和始末狀態軌道傾角均相差較大時，相對固定的擬合曲線導致形狀法很難收斂到可行解。

為了解決上述不足，本文在柱坐標系上建立了常值推力作用的航天器三維空間軌道機動動力學模型，在推力大小恒定的前提下，通過調整推力在空間中的俯仰角和方位角實現航天器在三維空間的軌道機動。本文采用正向設計法，提出了一種先軌道等待再軌道機動的優化設計策略，該策略能在軌道偏心率和始末狀態軌道傾角均相差較大的情況下實現航天器在三維空間的軌道機動任務。本文方法為三維空間軌道機動優化設計問題提供了新的思路。

1 軌道機動動力學模型

首先定義空間運動相關坐標系，如圖1 所示。定義EJ2000 坐標系OXIYIZI的原點O為地心，基準平面為歷元平赤道面，OXI軸由地心指向平春分點?，OZI軸與基準平面的法向平行，OYI軸的方向通過右手法則確定，OXIYIZI為靜系。定義參考坐標系ORXRYRZR的原點OR為航天器質心，ORZR軸與OZI軸平行，ORXR軸在矢量p和ORZR軸構成的平面內且與ORZR軸垂直，ORYR軸的方向通過右手法則確定，ORXRYRZR為動系。定義柱坐標系Orθz的原點為地心，極軸χ與OXI軸重合，以歷元平赤道面為基準面，r在基準面內且與極軸χ的夾角為θ，z與OZI軸平行。

圖1 空間運動坐標系之間關系Fig.1 Relationship between spatial motion coordinate systems

假設航天器的位置矢量p和速度矢量v在EJ2000 中的投影分別為pI=[pxI，pyI，pzI]T速度矢量為vI=[vxI，vyI，vzI]T=由圖1可知：

對式（1）求導可得：

根據位置矢量p的絕對導數與相對導數關系可得：

式中：S(ω)=為ω=[ω1，ω2，ω3]T的反對稱矩陣。dpdt表示在靜系中矢量p對時間t的絕對導數表示在動系中矢量p對時間t的相對導數，ω表示動系相對于靜系的角速度矢量。v和p在參考坐標系ORXRYRZR中的投影如圖2所示。

圖2 相關矢量在參考坐標系的投影Fig.2 Projections of the relevant vectors in the reference coordinate system

根據式（3）可得：

根據速度矢量v的絕對導數與相對導數關系可得：

式中：m為航天器的質量，p=‖p‖=為航天器質心到地心的距離，μ為地心引力常數，F為航天器所受控制力。F在參考坐標系ORXRYRZR中的投影如圖2 所示，其中φ為俯仰角，η為方位角。根據式（5）可得：

式中：F=‖F‖。航天器質量變化率為：

式中：ve表示航天器推進系統噴氣速度。根據式（4）、式（6）和式（7）可得到航天器軌道機動動力學模型：

式中：X(t)=(r，θ，z，vr，vt，vz)描述了航天器的狀態。只需要知道航天器在慣性系中描述的初始狀態和終端狀態，通過式（1）和式（2）將其轉換為柱坐標系中描述的初始狀態X(t0)和終端狀態X(tf)，利用優化算法對控制力F、俯仰角φ和為方位角η進行優化求解，便可實現航天器軌道機動優化設計。在工程實際中航天器的推力一般是恒定的，F和ve是常值，質量變化率也為常值。本文假設推力方向與航天器體軸平行，通過優化航天器飛行過程中的俯仰角φ和方位角η實現軌道機動優化設計。

2 參數無量綱化

式（8）描述的軌道機動動力學模型中各參數存在數量級差別，如果直接求解容易產生奇異，且影響計算效率。因此，需要對式（8）中各參數進行無量綱化處理。本文取航天器的初始質量m0為質量基本量綱，取地球半徑rE為長度基本量綱。力、速度和時間的基本量綱分別為和角度的基本量綱為1。式（8）中各參數經無量綱化處理之后的表達式為：

對式（9）～式（14）分別進行微分可得：

根據式（8）～式（21）可得到無量綱化的航天器軌道機動動力學模型：

3 優化設計策略

航天器軌道機動優化設計問題的實質是尋找一條從初始狀態到目標狀態的飛行軌跡，在滿足相關約束條件的同時使得航天器燃料消耗最低。傳統的軌道機動優化設計策略研究的是直接軌道機動問題，航天器從初始狀態到目標狀態的飛行過程中全程施加推力。本文提出一種先軌道等待再軌道機動的優化設計策略，航天器先在初始軌道上無動力飛行一段時間進行軌道等待，然后尋找一個最優時刻施加推力再進行軌道機動。與傳統的軌道機動優化策略相比，本文不僅需要尋找飛行過程中的最優控制量，還需要尋找最優等待時刻。本文優化設計的數學描述如下式所示：

圖3 優化設計策略Fig.3 Optimal design strategy

本文將多重直接打靶法（Multiple direct shooting method）和內點法相結合對優化問題進行具體求解。與直接打靶法相比，多重直接打靶法具有計算精度高且收斂性強的優點。其思路是將積分時間段劃分成若干個子時間區間，將各個子時間區間的待優化量進行離散化處理。已知初始狀態可以依次求出各子時間區間的狀態，中間每個子時間節點的狀態既是上一個子時間區間的末端狀態又是下一個子時間區間的起始狀態。這樣，將原先復雜的軌跡優化設計問題轉化為每個子時間區間的非線性規劃（Non-linear programming，NLP）問題，通過內點法求解得到最優飛行軌跡、最優俯仰角、最優方位角和

4 仿真算例

本文以航天器從近地小橢圓軌道異面機動轉移到遠地大橢圓軌道為例進行仿真分析，航天器進行軌道機動的始末端軌道根數如表1所示。航天器質量為3 000 kg，引力常數μ=3.986 × 1014m3s2，長度的基本量綱rU 為地球半徑。時間和速度的基本量綱分別為tU=806.785 6 s和vU=7 905.4 m s，質量的基本量綱為mU=3 000 kg，力的基本量綱為fU=29 396 N。無量綱化后的航天器推力和噴氣速度分別為0.1 fU 和6.202 2 vU。將表1 中機動軌道始末端軌道根數轉化為EJ2000 坐標系描述的始末端狀態，通過坐標轉換并進行無量綱化處理，得到表2所示的柱坐標系描述的無量綱化始末端狀態。

表1 機動軌道始末端軌道根數Table 1 Orbital elements at the beginning and end of the maneuvering orbit

表2 無量綱化始末端狀態Table 2 Dimensionless initial states and end states

本文設置了兩組仿真場景以驗證所提軌道機動優化設計策略，即：

1）先軌道等待再軌道機動

2）直接軌道機動

場景一：將軌道等待段和軌道機動段各劃分為10個子時間區間，每個子區間設置30個離散點。俯仰角和方位角的上界和下界分別設為π 和-π。在[ -π，π]范圍內，俯仰角和方位角的初值選取不影響多重直接打靶法求解的收斂性。本文中，各離散點的俯仰角和方位角初值均設為1。時刻和時刻的初值選取對多重直接打靶法求解的收斂性有影響，時刻的初值選取根據初始軌道周期來估計，時刻的初值選取根據終端軌道周期來估計。本文中，時刻和的初值分別設為5 和10，其上界分別設為7.884 2 和20.379 2，下界均設為0。仿真結果如圖4～圖7所示。

圖4 三維空間中的位置變化Fig.4 Position changes in three-dimensional space

圖4為航天器在三維空間中的位置變化，由圖4可知，航天器在初始軌道先飛行一段時間進行軌道等待，在t1時刻開始進行軌道機動，t2時刻到達目標位置。圖4 表明，在軌道偏心率和始末狀態軌道傾角均相差較大的情況下，航天器依然能實現三維空間的軌道機動任務。

圖5為航天器速度隨時間的變化關系，在[t0，t1)時間段不施加作用力，航天器進行無動力飛行。在[t1，t2]時間段施加常值推力，航天器進行軌道機動，總機動時間為5.788 6 tU，整個任務的總飛行時間為8.108 tU。

圖5 速度隨時間的變化Fig.5 Variation of velocity over time

圖6為航天器質量隨時間的變化關系，在[t0，t1)時間段航天器質量不變，在[t1，t2]時間段航天器質量勻速減小，在t2時刻航天器質量為0.906 6 mU。

圖6 質量隨時間的變化Fig.6 Variation of mass over time

圖7 為俯仰角和方位角隨時間的變化關系，在常值推力條件下，航天器通過控制俯仰角和方位角實現三維空間軌道機動。

圖7 俯仰角和方位角隨時間的變化Fig.7 Variations of pitch angle and azimuth angle over time

場景二：將軌道機動段劃分為10 個子時間區間，每個子時間區間設置30個離散點。俯仰角和方位角的設置與場景一相同。時刻的初值設為10，其上界和下界分別設為20 和0。仿真結果如圖8～圖11所示。將場景二與場景一進行對比可知，航天器在場景二中沒有軌道等待段，在初始時刻即進行軌道機動，tf時刻到達目標位置，實現了預定軌道機動任務。場景二的總機動時間為6.933 tU，機動時間比場景一多1.144 4 tU，任務總飛行時間比場景一少1.175 tU。到達目標位置時場景二中航天器質量為0.888 1 mU，場景二的燃料消耗比場景一多0.018 5 mU。場景二仿真結果表明，在完成相同的軌道機動任務情況下，本文所提軌道機動策略的燃料消耗更低。

圖8 三維空間的位置變化Fig.8 Variations of position in three-dimensional space

圖9 速度隨時間的變化Fig.9 Variation of velocity over time

圖10 質量隨時間的變化Fig.10 Variation of mass over time

圖11 俯仰角和方位角隨時間的變化Fig.11 Variations of pitch angle and azimuth angle over time

為了進一步分析航天器軌道機動過程中的燃料消耗，將本文方法與非連續推力的脈沖變軌方式進行對比。以脈沖施加時刻和相應的脈沖速度增量為待優化量，以總脈沖速度增量作為優化目標。各方向脈沖速度增量的上界和下界分別為0.05 vU 和-0.05 vU。根據場景二總燃料消耗折算成的總速度增量可知，最多施加10次脈沖便可完成預定的軌道機動任務，因此本文采用10次脈沖方式進行軌道機動。相鄰脈沖施加時刻之間的間隔上界為7.884 2，下界為0。本文中，各方向脈沖速度增量的初值均設為0.01 vU，相鄰脈沖施加時刻之間的間隔初值均設為1。軌道機動的始末端軌道參數和各物理量與連續推力情況相同，非連續推力的脈沖變軌仿真結果如圖12 所示。航天器的脈沖施加時刻和各方向速度增量如表3所示。

表3 脈沖時刻和速度增量Table 3 Pulse moments and velocity increments

圖12 三維空間的位置變化Fig.12 Variations of position in three-dimensional space

由圖12 可知，經過8 次脈沖變軌，航天器實現了預定的軌道機動任務，航天器總飛行時間為10.551 2 tU。將圖12 與圖4、圖8 進行對比可知，與連續推力方式的軌道機動不同，脈沖變軌方式先抬升軌道高度，在遠地點附近施加脈沖以改變軌道面。在完成軌道機動任務的前提下實現燃料優化目標。整個軌道機動任務過程中，航天器的總脈沖速度增量為0.491 3 vU，經換算，到達目標位置時航天器質量為0.923 8 mU。非連續推力的脈沖變軌方式燃料消耗比場景一低0.017 2 mU，這是因為連續推力變軌過程中存在無效推力而造成額外質量損失，脈沖變軌方式是理想化的變軌方式。

為了驗證本文方法的在軌道機動任務中的普適性，有必要對航天器在平面內的軌道機動進行仿真分析。平面內軌道機動的始末端軌道根數如表4所示，表5 為無量綱化處理后的柱坐標系描述的始末端狀態。各仿真參數設置與場景一相同，仿真結果如圖13～圖16所示。

表4 機動軌道始末端軌道根數Table 4 Orbital elements at the beginning and end of the maneuvering orbit

表5 無量綱化始末端狀態Table 5 The dimensionless initial states and end states

圖13 三維空間中的位置變化Fig.13 Variations of position in three-dimensional space

圖14 速度隨時間的變化Fig.14 Variation of velocity over time

圖15 質量隨時間的變化Fig.15 Variation of mass over time

圖16 俯仰角和方位角隨時間的變化Fig.16 Variations of pitch angle and azimuth angle over time

由圖13可知，航天器在初始軌道先飛行一段時間進行軌道等待，在t1時刻開始進行軌道機動，t2時刻到達目標位置，實現了平面內的軌道機動任務。航天器無動力飛行時間為3.16 tU，總機動時間為5.917 tU，在t2時刻航天器質量為0.906 6 mU。

5 結論

本文研究了常值推力作用的航天器軌道機動優化設計問題。建立了無量綱化的軌道機動動力學模型并提出了先軌道等待再軌道機動的優化設計策略。通過多重直接打靶法將復雜的軌跡優化問題轉化為非線性規劃問題并采用內點法進行求解。仿真結果表明，本文所提方法在軌道機動任務中的具有普適性，不僅適用于平面軌道機動任務，同樣適用于空間軌道機動任務。與傳統方法相比，本文所提的軌道機動優化設計策略燃料消耗更低。