考慮J2項攝動的Lambert問題狀態(tài)空間攝動解析近似法

李兆亭 張洪波 鄭 偉

國防科技大學(xué)空天科學(xué)學(xué)院，長沙 410073

Lambert問題是航天動力學(xué)中的經(jīng)典兩點邊值問題，在航天器軌道設(shè)計、軌道確定、軌道機(jī)動與制導(dǎo)、彈道導(dǎo)彈的瞄準(zhǔn)射擊和顯式制導(dǎo)等方面有廣泛的應(yīng)用。

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

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

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

1 Lambert問題基本模型

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

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

(1)

根據(jù)x可求得轉(zhuǎn)移軌道的半通徑p：

(2)

初始點處的轉(zhuǎn)移軌道速度：

(3)

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

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

Δv=v1-v0

(4)

由此可得二體條件下的初始點處脈沖速度和標(biāo)準(zhǔn)的轉(zhuǎn)移軌道。

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

2 狀態(tài)空間攝動法

建立當(dāng)?shù)厮阶鴺?biāo)系O-rβz，坐標(biāo)系的原點位于航天器的質(zhì)心，r軸沿地心位置矢量r的方向；β軸在軌道面內(nèi)垂直于r軸，沿飛行方向為正；z軸由右手法則確定，沿動量矩h的反方向。該坐標(biāo)系可記為LVLH坐標(biāo)系，如圖1所示。

圖1 LVLH坐標(biāo)系

(5)

(6)

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

當(dāng)航天器在LVLH坐標(biāo)系下的三軸上有攝動加速度時，實際運(yùn)動參數(shù)將偏離預(yù)定二體軌道。

取航天器運(yùn)動的狀態(tài)向量為：

(7)

其二體運(yùn)動狀態(tài)方程(6)可以寫成如下形式：

(8)

設(shè)該方程的解為X0=X0(β)。考慮攝動影響后，方程變?yōu)椋?/p>

(9)

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

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

(10)

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

(11)

本問題中，轉(zhuǎn)移軌道的時間是一定的，故要研究2種運(yùn)動的等時偏差。將其運(yùn)動偏差狀態(tài)矢量取為：

(12)

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

(13)

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

根據(jù)線性系統(tǒng)的性質(zhì)，系統(tǒng)的狀態(tài)響應(yīng)可以表示為零輸入響應(yīng)Δt,u0X和零狀態(tài)響應(yīng)Δt, x0X之和，即：

(14)

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

(15)

3 基于狀態(tài)空間攝動法的初始速度修正

3.1 修正速度計算

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

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

(16)

同時在任務(wù)開始時刻，初始位置也是給定的，無法瞬時調(diào)整，即3個零輸入響應(yīng)的位置參量為0：

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

(17)

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

(18)

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

(19)

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

3.2 零狀態(tài)響應(yīng)偏差量計算

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

(20)

等角攝動偏差量與等時偏差攝動量滿足以下關(guān)系式：

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

(21)

建立的線性時變方程為：

(22)

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

若以起始點K作為地心航程角的起算點，則等角攝動零狀態(tài)響應(yīng)可表示為：

(23)

求得等角攝動的零狀態(tài)響應(yīng)之后，根據(jù)線性系統(tǒng)的性質(zhì)可計算得到等時零狀態(tài)響應(yīng)為：

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

(24)

下面給出等角攝動零狀態(tài)響應(yīng)的解析解。

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

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

(25)

側(cè)向偏差解析解為：

(26)

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

(27)

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

(28)

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

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

(29)

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

圖3 超幾何級數(shù)的近似效果

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

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

4 仿真分析

4.1 仿真算例

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

表1 基本Lambert問題參數(shù)

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

表2 轉(zhuǎn)移軌道起點處的軌道根數(shù)比較

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

加入修正速度后，考慮J2項進(jìn)行軌道積分得到終點的位置速度，結(jié)果如表3所示，其大小為19m?？梢奐2項的影響得到了很好的修正。

表3 位置偏差

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

表4 狀態(tài)空間攝動方法與傳統(tǒng)迭代方法的比較結(jié)果

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

為深入分析該方法的性能，本小節(jié)將通過改變轉(zhuǎn)移軌道的軌道根數(shù)來觀察狀態(tài)空間攝動法的修正效果。設(shè)標(biāo)準(zhǔn)轉(zhuǎn)移軌道的根數(shù)如表5所示。

表5 轉(zhuǎn)移軌道的標(biāo)準(zhǔn)軌道根數(shù)

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

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

其次，改變轉(zhuǎn)移軌道的軌道傾角，結(jié)果如圖5所示?？梢姡壍纼A角在[0°，180°]的范圍內(nèi)，其偏差總體上為兩端大中間小，在90°時有極小值，位置偏差均優(yōu)于10m。同樣，偏差隨偏心率增大而增大。

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

改變起點真近點角的結(jié)果如圖6所示?？梢?，真近點角在[-180°，180°]的范圍內(nèi)，其偏差總體上為兩端小中間大，在-45°～0°之間存在極大值，位置偏差均優(yōu)于10m。

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

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

圖7 位置偏差隨轉(zhuǎn)移時間的變化

5 結(jié)論

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

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