考慮J2項攝動的Lambert問題狀態空間攝動解析近似法

李兆亭 張洪波 鄭 偉

國防科技大學空天科學學院，長沙 410073

Lambert問題是航天動力學中的經典兩點邊值問題，在航天器軌道設計、軌道確定、軌道機動與制導、彈道導彈的瞄準射擊和顯式制導等方面有廣泛的應用。

在二體條件下，Lambert問題已經有許多經典的解法[1-3]。在實際任務中，如航天器軌道轉移、遠程彈道導彈制導等，地球非球形二階帶諧項(J2項)的影響不可忽略，需要對其進行修正。目前已經有多種修正方法[4-6]，其中最經典的方法是利用二體狀態轉移矩陣和微分方程不斷迭代修正初始速度，以滿足終端的位置精度要求。如黨朝輝[7]等提出了一種基于狀態傳遞矩陣解析解的J2攝動Lambert問題求解方法，采用經典的二體Lambert算法求解轉移軌道初始速度，然后采用循環迭代的計算流程將其修正到考慮J2攝動的條件下。該方法計算效率高、速度快。呼衛軍[8]等圍繞參考攔截面內與法向方向攝動展開項，結合理想引力場解算公式，改進了傳統的需求速度求解方法，所需計算量小、計算精度高。孫瑜[9]等利用平根數法解攝動方程，提出了一種考慮J2攝動并包含短周期項的解析法，并將其用于彈道預報。魏倩[10]等提出了一種基于神經網絡的改進制導方法，該方法利用BP神經網絡設計了一種J2項攝動的虛擬目標點預測模型并補償攝動偏差，其結果精度高、計算規模小。Zhang[6]等基于打靶法的原理提出了一種新的偏移迭代法，其收斂速度快、與硬件結合好。Marcus[11]等提出了一種J2項攝動下使用高斯變分方程計算軌道運行范圍的方法，同時解決了J2攝動下的自由時間最小脈沖軌道轉移問題。Robyn[12]等提出了一種使用Kustaanheimo-Stiefel變換和改進的切比雪夫-皮卡德迭代解決兩點邊值問題和初值問題的新方法，用解析方法解決了攝動條件下的兩點邊值問題。Roberto[13]等針對多圈、擾動條件下的Lambert問題，提出了2種基于高階泰勒展開式的求解器，并在此基礎上求出了考慮J2項攝動問題的解析解。Woollands[14]等提出了使用特殊解法和修改的切比雪夫-皮卡爾迭代法來解決多圈攝動Lambert問題的方法。該特殊解法使用參考軌跡和一組特定解，結合修改的切比雪夫-皮卡爾迭代法，提高了算法的效率，保持了計算精度。

任萱[15]、鄭偉[16]等提出的狀態空間攝動法是一種將非線性方程線性化從而得到近似解析解的方法。其基本思路為：基于一條標準彈道將考慮攝動的非線性方程線性化，得到以位置速度偏差為狀態變量的線性系統，求解其狀態轉移矩陣即可得到通解，進而得到近似解析解。它在彈道導彈主動段、自由端彈道攝動分析等方面有良好的效果[16]。

基于狀態空間攝動法，提出一種考慮J2項攝動的Lambert問題解法。該方法利用建立的線性化攝動運動方程，令線性系統的零輸入響應和零狀態響應在目標點處相互抵消，從而求得初始速度修正量的解析表達式，最后通過數值仿真結果驗證了方法的有效性。

1 Lambert問題基本模型

Lambert問題是給定初始點、終端點的位置矢徑以及轉移時間，求滿足位置和時間要求的轉移軌道的問題。常見的解決Lambert問題的方法有：高斯迭代法、P迭代法、普適變量法及Battin法等[1]。使用Battin法來求解Lambert問題[2-5]：

Battin法以x為自變量，通過迭代高斯第一和第二方程得到Lambert問題的解，x的表達式為：

(1)

根據x可求得轉移軌道的半通徑p：

(2)

初始點處的轉移軌道速度：

(3)

式(1)～(3)中,r1,r2為地心距；Δβ為地心轉移角；p為半通徑；E,x,y,m,s,λ為中間量[2]。

設v0為變軌前的速度，則初始點處的脈沖速度為：

Δv=v1-v0

(4)

由此可得二體條件下的初始點處脈沖速度和標準的轉移軌道。

考慮攝動的條件下，初始脈沖速度以式(4)所求為準，所得轉移軌道逐漸偏離標準軌道，終端點位置會有較大的偏差。本文采用狀態空間攝動法對此位置偏差進行修正，得到了修正后的脈沖速度。

2 狀態空間攝動法

建立當地水平坐標系O-rβz，坐標系的原點位于航天器的質心，r軸沿地心位置矢量r的方向；β軸在軌道面內垂直于r軸，沿飛行方向為正；z軸由右手法則確定，沿動量矩h的反方向。該坐標系可記為LVLH坐標系，如圖1所示。

圖1 LVLH坐標系

(5)

(6)

將航天器的二體軌道運動方程投影到LVLH坐標系中，得到的結果為式(5)。將式(5)的自變量由t變換為β，可得式(6)。

當航天器在LVLH坐標系下的三軸上有攝動加速度時，實際運動參數將偏離預定二體軌道。

取航天器運動的狀態向量為：

(7)

其二體運動狀態方程(6)可以寫成如下形式：

(8)

設該方程的解為X0=X0(β)。考慮攝動影響后，方程變為：

(9)

考慮到攝動量U(X,β)是小量, 上述方程的解可以寫成如下形式：

X=X0(β)+ΔX(β)

(10)

將式(10)代入到考慮攝動的狀態方程式(9)中，將其線性化展開并略去二階以上的小量，可得到運動偏差量的線性化方程：

(11)

本問題中，轉移軌道的時間是一定的，故要研究2種運動的等時偏差。將其運動偏差狀態矢量取為：

(12)

綜合式(11)和(12)，建立的線性時變狀態方程[15-16]如下：

(13)

式中，矩陣D，C,U的表達式和式(13)的狀態轉移矩陣Φt(β1,β0)詳見文獻[15-16]。

根據線性系統的性質，系統的狀態響應可以表示為零輸入響應Δt,u0X和零狀態響應Δt, x0X之和，即：

(14)

式中，ΔtX(β1)為考慮攝動條件下的終端點處的狀態偏差量，Δt, u0X(β1)表示由系統初始狀態轉移得到的零輸入狀態偏差量，Δt, x0X(β1)表示在攝動影響下的零狀態運動偏差量，ΔtX(β0)為初始點處的狀態偏差量，U(ξ)為攝動項。上式展開容易得到：

(15)

3 基于狀態空間攝動法的初始速度修正

3.1 修正速度計算

根據Lambert問題的定義，終端點的位置是給定的，因此要求有：

Δtr(β1)=0, Δtβ(β1)=0, Δtz(β1)=0

(16)

同時在任務開始時刻，初始位置也是給定的，無法瞬時調整，即3個零輸入響應的位置參量為0：

Δt,u0r(β0)=0, Δt,u0β(β0)=0, Δt,u0z(β0)=0

(17)

將式(16)和(17)代入到式(15)中可得：

(18)

考慮上述公式中Φt矩陣的第2、4和6行所在的等式，只需要知道零狀態響應下的3個位置偏量，即可求出起始點處的修正速度：

(19)

式中,Φt已知，故問題轉換為計算Δt,x0r，Δt,x0β和Δt,x0z三個終端點處的零狀態偏差量。

3.2 零狀態響應偏差量計算

由于直接計算等時攝動零狀態響應方程時，不便于確定終端范圍角β1。因此，可以先計算等角攝動的零狀態響應，然后通過變換矩陣將其轉換為等時攝動零狀態響應[15]。等角攝動的運動偏差狀態量取為：

(20)

等角攝動偏差量與等時偏差攝動量滿足以下關系式：

ΔrX(β)=D(β)·ΔβX(β)

(21)

建立的線性時變方程為：

(22)

矩陣C，D，J2項引力攝動加速度U的表達式和方程(22)的狀態轉移矩陣Φβ(β1,β0)有解析解，結果見文獻[15-16]。

若以起始點K作為地心航程角的起算點，則等角攝動零狀態響應可表示為：

(23)

求得等角攝動的零狀態響應之后，根據線性系統的性質可計算得到等時零狀態響應為：

Δr, x0X(β1)=D(β1)·Δβ, x0X(β1)

(24)

下面給出等角攝動零狀態響應的解析解。

3.3 等角攝動偏差量的解析解

根據文獻[17]所提的極點變換方法，重新定義航程角，并將J2項攝動加速度分解為縱向和側向加速度，分別求出其等角條件下的攝動解析解。其中，縱向偏差解析解為：

(25)

側向偏差解析解為：

(26)

同時，由于偏差解析解的積分中涉及使用了偏近點角E，其與真近點角關系為：

(27)

式(25)和(26)中需要對真近點角進行積分，故使用了E對于f的超幾何級數對式(27)進行近似，取前5階的結果如式(28)所示。

(28)

對于式(27)，當f∈[-π,+2π]時，為了保證積分的連續性，使得Ek≥E0(終點大于起點的偏近點角)，取值如式(29)所示。

在f∈[-π,+2π]，e∈[0,0.5]時，式(29)得到的實際偏近點角如圖2所示，式(28)與(29)得到的偏近點角之差如圖3所示。

(29)

圖2 偏近點角隨真近點角和偏心率的變化

圖3 超幾何級數的近似效果

從圖2和3中可以看到，對于f∈[-π,+2π],e∈[0,0.5]范圍內的數據，相對于達到幾百度的實際偏近點角值的偏差最多不超過0.01°，式(28)的近似效果較好。

根據式(25)～(28)即可求出該解析解，其中文獻[18]給出了等角攝動零狀態響應解析解的4個值，分別為Δβvr，Δβr，Δβvβ和Δβz。關于時間的Δβt和關于法向速度的Δβvz也可通過推導得到解析解。

4 仿真分析

4.1 仿真算例

仿真分析中始末點的條件和飛行時間如表1所示。

表1 基本Lambert問題參數

根據Battin的算法可以得到二體條件下的轉移軌道，軌道根數如表2所示。以二體軌道的初始速度開始，考慮J2項進行軌道積分可以得到終點的位置速度，終端位置偏差如表3所示。可見，位置偏差為18.065km。

表2 轉移軌道起點處的軌道根數比較

根據本文所提的修正方法，可求出考慮J2項攝動后初始速度的修正量，在LVLH坐標系中其大小為：

加入修正速度后，考慮J2項進行軌道積分得到終點的位置速度，結果如表3所示，其大小為19m。可見J2項的影響得到了很好的修正。

表3 位置偏差

同時，該方法與傳統的二體轉移矩陣迭代計算的比較結果如表4所示。其中，表中的Δvx，Δvy和Δvz為起點處修正后的速度增量，2種方法差別極小，小于0.01m/s。Δr為終端時刻的位置偏差，可見狀態空間攝動法得到的終端位置精度較差，這與其方法本身的近似線性化有關。T0為2者在相同條件下的仿真計算時間，相比狀態空間攝動法，不需要迭代積分過程，因而計算速度更快。

表4 狀態空間攝動方法與傳統迭代方法的比較結果

4.2 不同始末條件下的修正效果比較

為深入分析該方法的性能，本小節將通過改變轉移軌道的軌道根數來觀察狀態空間攝動法的修正效果。設標準轉移軌道的根數如表5所示。

表5 轉移軌道的標準軌道根數

首先，改變轉移軌道的半長軸和偏心率，結果如圖4所示。可見，修正后的終端位置偏差都優于20m。同時，軌道半長軸影響了修正后的位置偏差，而且軌道高度越高，其修正偏差越小,與J2項攝動隨高度變化的影響是一致的。同時，偏心率越大，其偏差越大;

圖4 位置偏差隨軌道長半軸的變化

其次，改變轉移軌道的軌道傾角，結果如圖5所示。可見，軌道傾角在[0°，180°]的范圍內，其偏差總體上為兩端大中間小，在90°時有極小值，位置偏差均優于10m。同樣，偏差隨偏心率增大而增大。

圖5 位置偏差隨軌道傾角的變化

改變起點真近點角的結果如圖6所示。可見，真近點角在[-180°，180°]的范圍內，其偏差總體上為兩端小中間大，在-45°～0°之間存在極大值，位置偏差均優于10m。

圖6 位置偏差隨初始真近點角的變化

改變轉移時間的結果如圖7所示。T為轉移軌道的周期。終端位置偏差隨轉移時間的增長而增長，這主要是由狀態空間攝動法的線性化誤差引起的。位置偏差都優于300m，修正效果較好。在位置偏差最大時，地心轉移角約為227.88°。

圖7 位置偏差隨轉移時間的變化

5 結論

研究了一種基于狀態空間攝動法的考慮J2項下的Lambert問題的修正方法，并對修正效果的影響因素進行了分析。該方法主要利用建立好的線性化攝動運動方程，令線性系統的零輸入響應和零狀態響應在目標點處相互抵消，從而得到修正量的解析表達。

仿真結果表明，該方法對于考慮J2項下的Lambert問題有很好的修正效果，位置精度得到了很大的提高，但同時其精度隨轉移時間的增加而變低。該方法可用于軌道機動與制導、彈道導彈的瞄準射擊等方面。