高精度空間非合作式相對軌道狀態(tài)預(yù)報*

朱正龍，閆　野，楊躍能

(國防科技大學 航天科學與工程學院， 湖南 長沙　410073)

?

高精度空間非合作式相對軌道狀態(tài)預(yù)報*

朱正龍，閆野，楊躍能

(國防科技大學 航天科學與工程學院， 湖南 長沙410073)

摘要：為了解決接近空間非合作目標過程中的相對導航問題，提出一種基于相對軌道動力學方程和二階龍格庫塔積分公式的高精度空間非合作目標相對軌道預(yù)報模型?？紤]地球J2攝動加速度，將目標軌道方程在參考軌道附近展開，保留至引力加速度差的二階展開項，并進一步推導考慮J2攝動的參考軌道角速度和角加速度，建立相對軌道動力學微分方程；在給定相對軌道初值的情況下，采用二階龍格庫塔積分公式進行相對軌道預(yù)報。該模型采用數(shù)值積分方法，對相對軌道動力學模型形式?jīng)]有限制，通用性好，適用范圍廣；選擇低階龍格庫塔公式積分預(yù)報，既減少了計算量又保證了計算精度。設(shè)置高軌和低軌兩種相對運動仿真場景，仿真結(jié)果表明了模型的通用性和精確性。

關(guān)鍵詞：相對軌道動力學；J2攝動；相對軌道預(yù)報

空間相對軌道狀態(tài)預(yù)報技術(shù)是完成空間接近操作任務(wù)的關(guān)鍵技術(shù)之一。針對空間非合作目標的相對軌道預(yù)報難度更大，不僅需要追蹤航天器對非合作目標進行跟蹤測量，還需要不依賴目標軌道參數(shù)的相對軌道預(yù)報模型提供理論支撐。

經(jīng)典相對軌道預(yù)報模型，如CW模型和TH模型[1-2]，廣泛應(yīng)用于空間交會對接、編隊飛行以及其他在軌操作接近任務(wù)中，其優(yōu)點是可獲得解析解，使用方便，但是該模型忽略了地球非球形引力攝動，并對相對運動方程進行了線性化處理，適用范圍有限，預(yù)報精度也受到限制，難以滿足高精度空間操作任務(wù)的需求。Schweighart和Sedwick[3]考慮J2攝動的影響，并引入平均參考軌道的概念，建立了線性化的相對軌道模型，適用于近圓軌道的空間編隊任務(wù)。該模型雖然也可以獲得解析解，但是參考軌道的平均化、相對運動方程的線性化處理以及近圓假設(shè)均限制了模型的適用范圍。Pluym和Damaren[4]建立了考慮中心引力3階、J2攝動引力2階的相對運動模型，該模型僅給出微分方程，并且采用了平均軌道角速度。Chen和Jing[5]基于拉格朗日動力學理論，建立了考慮J2攝動和大氣阻力攝動的精確相對運動微分方程，并給出了相應(yīng)的精確軌道角速度矢量和軌道角加速度矢量。

現(xiàn)有的相對運動模型存在的問題是適用范圍受限、通用性差或者難以求解，這都不利于實際工程應(yīng)用。

1相對軌道動力學

1.1一般相對動力學方程

為描述問題的方便，首先定義軌道坐標系：原點o位于航天器的質(zhì)心，ox軸由地心指向航天器的質(zhì)心，oy在軌道面內(nèi)垂直于ox軸并指向運動方向， ox，oy和oz構(gòu)成右手直角坐標系。

設(shè)追蹤航天器和目標航天器的位置矢量分別為Rc和Rt，追蹤航天器的軌道角速度矢量為ωo。定義兩個航天器的相對位置矢量為ρ=Rt-Rc，根據(jù)動坐標系下的矢量微分法則，在追蹤航天器的軌道坐標系中：

(1)

(2)

式中，Δgc表示中心引力加速度差，Δgp為攝動加速度差之和。聯(lián)立式(1)和式(2)可得追蹤器軌道坐標系下的目標相對運動方程為：

(3)

1.2精確二體相對運動模型

在式(3)中，忽略攝動加速度和控制加速度，可得到精確的二體相對運動方程：

(4)

(5)

在二體運動條件下，追蹤器軌道角速度矢量和軌道角加速度矢量可表示為[1]：

(6)

式中，θ為真近點角。目標器和追蹤器的中心引力加速度在追蹤器軌道坐標系中的投影分別為：

(7)

(8)式中，μ是地球引力常數(shù)，rt為目標器的地心距，rc為追蹤器的地心距。中心引力加速度差的表達式為：

(9)

(10)

式(10)即為二體條件下目標相對追蹤器的精確運動方程，為時變系數(shù)的非線性微分方程組。

1.3精確J2攝動相對運動模型

文獻[5]在推導相對軌道微分方程時，考慮了J2攝動和大氣阻力攝動，本文僅考慮J2項地球扁率攝動，忽略其他攝動加速度和控制加速度，建立相對運動模型：

此時，目標器與追蹤器的軌道均受到J2攝動力的影響。根據(jù)文獻[6]，攝動加速度在軌道坐標系中可表示為：

(12)

式中，r表示地心距(即航天器質(zhì)心與地心之間的距離)，i為軌道傾角，u為緯度幅角，RE是地球半徑，J2是地球引力位系數(shù)。

在追蹤器軌道坐標系中，地球扁率攝動加速度差的計算表達式為：

(13)

軌道角速度矢量ωo=[ωx,ωy,ωz]T的計算過程及表達式為：

(14)

式中，Ω為升交點赤經(jīng)，Mj(·)(j=1,2,3)表示繞坐標軸的初等轉(zhuǎn)換矩陣。進一步整理式(14)可得：

(15)

軌道根數(shù)的導數(shù)可根據(jù)高斯攝動方程獲得(參見文獻[7])，將式(12)代入高斯攝動方程可得[5]：

(16)

(17)

(18)

(19)

(20)

(21)

將式(9)、式(13)、式(20)和式(21)代入式(11)即得到目標器相對追蹤器的軌道運動模型。

2相對運動模型的展開

從式(13)可以看出，在計算J2攝動加速度差時，還需要獲取目標航天器的軌道參數(shù)，而對非合作空間目標，其軌道參數(shù)常常無法獲取。為此考慮將目標航天器軌道在參考軌道航天器附近進行泰勒級數(shù)展開。

2.1中心引力加速度的展開

根據(jù)式(7)和式(8)，將gc,t在gc,c附近級數(shù)展開可得：

(22)

式中，ei為三維向量空間中的標準基向量，HOT1表示目標引力加速度gc,t(ρ)關(guān)于ρ的高階展開項。在式(22)中，作如下符號定義：

P1(rc)

(23)

Q1i(rc)

(24)

經(jīng)計算推導，可得：

(25)

(26)

(27)

(28)

根據(jù)式(22)至式(28)，經(jīng)推導計算，中心引力加速度差的一階項和二階項表達式分別為：

(29)

(30)

2.2J2攝動加速度的展開

根據(jù)式(12)，目標J2攝動加速度gJ2,t在追蹤器軌道坐標系中的計算式為：

(31)

為方便計算，將gJ2,t轉(zhuǎn)化為顯含ρ的函數(shù)表達式。根據(jù)文獻[1]，目標器在慣性坐標系下的J2項攝動加速度可表示為：

(32)

將Rt=[X,Y,Z]T投影到追蹤器軌道坐標系，其關(guān)系式為：

(33)

(34)

(35)

將式(32)至式(34)代入式(35)，可將gJ2,t表示為顯含相對位置矢量ρ的函數(shù)gJ2,t(ρ)，將gJ2,t(ρ)在gJ2,c附近級數(shù)展開可得：

(36)

式中，HOT2表示目標引力加速度gJ2,t(ρ)關(guān)于ρ的高階展開項。作如下符號定義：

P2(rc,ic,uc)

(37)

Q2i(rc,ic,uc)

(38)

記sinφ=sφ,cosφ=cφ，φ為任意的角變量，經(jīng)過推導計算可得[4]：

(39)

(40)

(41)

(42)

式(39)～(42)省略了ic和uc的下標。根據(jù)式(36)，保留至J2攝動加速度差的二階近似表達式為：

將式(43)中的一階項和二階項分別記為ΔgJ2,1st和ΔgJ2,2nd，即：

ΔgJ2,1st=P2(rc,uc,ic)ρ

(44)

(45)

取中心引力和J2攝動加速度差的一階展開項，根據(jù)式(11)可得線性化的相對運動方程為：

(46)

(47)

式中，時變矩陣A(t)的表達式為：

(48)

方程式(11)給出了考慮J2地球扁率攝動的精確相對運動方程，其中J2攝動加速度差的計算包含了目標航天器的軌道根數(shù)，這對于非合作目標的相對軌道計算很不方便。為此，考慮取全部中心引力加速度差和J2攝動加速度差的前兩階展開項，得到高精度的非線性相對運動方程為：

(49)

式中，Δgc和ΔgJ2,1,2分別由式(9)和式(43)給出。

3相對軌道預(yù)報模型

本節(jié)在相對軌道運動微分方程式(49)的基礎(chǔ)上，研究快速、高精度的空間非合作目標相對軌道狀態(tài)預(yù)報方法。

在進行相對軌道預(yù)報時，一般根據(jù)當前時刻的狀態(tài)x(k)及追蹤器軌道參數(shù)σc計算出經(jīng)h=Δt時間之后的狀態(tài)x(k+1)，可表示為如式(50)所示的非線性函數(shù)的形式：

xk+1=F(tk,xk)

(50)

為實現(xiàn)相對狀態(tài)的時間轉(zhuǎn)移，現(xiàn)在的任務(wù)就是構(gòu)造一個計算速度快、計算精度損失小的非線性矢量函數(shù)F:xkxk+1。借助數(shù)值積分的方法，采用低階的Runge-Kutta積分方法實現(xiàn)相對狀態(tài)的轉(zhuǎn)移。設(shè)微分方程形式的相對軌道模型可表示為：

(51)

式中，σc表示追蹤器的軌道根數(shù)，在此假設(shè)軌道根數(shù)σc,k可由追蹤器上的導航設(shè)備精確提供，并且在單個狀態(tài)轉(zhuǎn)移周期Δt內(nèi)保持不變。

3.1線性化的相對軌道預(yù)報模型

線性系統(tǒng)理論已經(jīng)十分成熟，而且在使用中十分方便。若線性化的相對軌道運動方程為：

(52)

則根據(jù)常微分方程的知識，求解上述連續(xù)時變系統(tǒng)，可得：

x(t)=Φ(t,t0)x(t0)

(53)

式中，Φ(t,t0)是狀態(tài)轉(zhuǎn)移矩陣，其計算公式為：

Φ(t,t0)=exp{A(t)(t-t0)}

(54)

對于矩陣的指數(shù)運算，其定義為：

(55)

從而，當狀態(tài)轉(zhuǎn)移時間Δt=t-t0較小時，可得到狀態(tài)轉(zhuǎn)移矩陣的表達式為：

(56)

考慮地球中心引力加速度差和J2攝動加速度差的一階展開模型由式(47)給出，相應(yīng)的系統(tǒng)矩陣A(t)由式(48)計算，階數(shù)根據(jù)需要選擇。

3.2非線性相對軌道預(yù)報模型

對于高精度的狀態(tài)預(yù)報問題，狀態(tài)方程線性化處理損失了模型的精度，為此須對非線性相對軌道運動方程進行求解。相對軌道運動方程由式(49)給出，記x=[ρT,vT]T，可將式(49)寫成式(51)的形式，并且：

f(x,t;σc)=

(57)

基于二階Runge-Kutta公式的相對狀態(tài)轉(zhuǎn)移公式為：

(58)

式中，h為積分步長。二階Runge-Kutta公式為二階O(h2)精度的積分公式，具體的狀態(tài)轉(zhuǎn)移精度由狀態(tài)轉(zhuǎn)移周期Δt決定。

4仿真分析

本節(jié)從積分方法、積分步長以及相對動力學方程3個方面分析所建立的相對軌道預(yù)報模型的精度。在仿真分析時，采用STK-HPOP模塊生成目標航天器和追蹤航天器的軌道數(shù)據(jù)，并轉(zhuǎn)換為相對軌道數(shù)據(jù)，作為真實軌道數(shù)據(jù)用于對比分析模型精度。

在高軌和低軌分別設(shè)計一個目標器與追蹤器先接近后遠離的相對運動場景，最小相對距離為20km，仿真時間為以目標器、追蹤器最近距離為中點的前后共600s的時間。

低軌仿真場景：目標器與追蹤器軌道分別運行于800km高度的太陽同步軌道上，兩軌道面的夾角為7.69°，目標器和追蹤器先后過軌道交點，最小相對距離約為20km。

高軌仿真場景：目標器運行于地球同步軌道上，追蹤器運行于周期為12小時、遠地點高度低于目標軌道20km的橢圓軌道上。

4.1積分方法對預(yù)報精度的影響

采用基于二階Runge-Kutta積分公式的相對軌道預(yù)報模型，分別在低軌和高軌兩個不同場景條件下，仿真分析考慮不同展開項引力加速度相對運動模型的預(yù)報精度。

圖1給出了分別基于Euler方法和二階至四階Runge-Kutta方法(RK2，RK3，RK4)進行軌道預(yù)報的誤差曲線，仿真計算步長為0.1s。仿真結(jié)果表明：采用二階Runge-Kutta方法既保證了計算精度，也提高了計算效率。

(a) 低軌算例(a) Low earth orbit example

(b) 高軌算例(b) High earth orbit example圖1　采用不同積分方法時的相對軌道預(yù)報誤差曲線Fig.1　Relative orbit propagation error withdifferent integral methods

4.2積分步長對預(yù)報精度的影響

將積分步長分別取為0.1s，1s，5s和10s四種情況分別進行軌道預(yù)報，經(jīng)過600s后，預(yù)報誤差如表1和表2所示。

表1和表2的結(jié)果表明：若采用二階以上(含二階)積分方法，模型中計算步長在0.1～10s之間時的誤差發(fā)散情況相當；而采用Euler方法時，誤差近似呈線性發(fā)散，計算精度很低。

表1　采用不同積分步長時的低軌相對軌道預(yù)報誤差

表2　采用不同積分步長時的高軌相對軌道預(yù)報誤差

4.3相對動力學方程對預(yù)報精度的影響

當相對軌道動力學方程保留不同階數(shù)的地球引力加速度差時，得到的相對軌道預(yù)報模型精度不同。圖2(a)和圖2(b)分別為低軌和高軌場景下的相對軌道預(yù)報誤差曲線。

(a) 低軌算例(a) Low earth orbit example

(b) 高軌算例(b) High earth orbit example圖2　采用不同動力學方程時的相對軌道預(yù)報誤差曲線Fig.2　Relative orbit propagation errors with differentequations of relative orbit dynamics

仿真結(jié)果表明：在600s的相對運動過程中，二體相對模型TB及其一階展開式TB(1)、二階展開式TB(1,2)在低軌場景下預(yù)報誤差均在50m以上，而在高軌的預(yù)報誤差小于2m；考慮J2攝動及其一階、二階展開項的模型，即TB+J2，TB+J2(1)和TB+J2(1,2)，在低軌的預(yù)報誤差均在10m以內(nèi)，而在高軌的預(yù)報誤差在1m以內(nèi)；此外兩個線性模型，即TB(1)和TB(1)+J2(1)在高軌的預(yù)報誤差在2m以內(nèi)，在低軌的預(yù)報誤差均超過1000m。

綜上可知：在高軌，采用二體模型的一階展開式(即著名的Lawden方程)建立預(yù)報模型即可獲得較高的預(yù)報精度，預(yù)報誤差在2m以內(nèi)；在低軌，除考慮中心引力差的一階展開項，還必須考慮中心引力的二階展開項以及J2攝動加速度差的一階展開項，才能獲得預(yù)報誤差小于10m的模型，若進一步考慮J2攝動加速度差的二階展開項，600s內(nèi)的預(yù)報誤差可縮小到2m以內(nèi)。

J2攝動是最主要的軌道攝動因素之一，考慮到J2攝動加速度的大小近似與航天器地心距的5次方成反比，而高軌航天器的地心距約為低軌的6倍，由此可知高軌J2攝動加速度的影響很小，而低軌影響明顯。此外，由于相對距離遠，攝動加速度差的高階展開項不再是小量，如果忽略會帶來較大的預(yù)報誤差。

5結(jié)論

空間非合作目標相對軌道預(yù)報是實現(xiàn)空間操作的一項重要技術(shù)環(huán)節(jié)?？紤]主要軌道攝動因素，并基于二階Runge-Kutta方法，提出了一種相對軌道預(yù)報模型，綜合考慮計算精度和計算效率，該模型通用性好，計算速度快，預(yù)報精度高，可用于高速相對運動空間非合作目標相對軌道預(yù)報，且在算例給定的時間內(nèi)，其預(yù)報誤差在米量級。所提方法能夠?qū)鉀Q未來更高精度的空間操作任務(wù)提供經(jīng)驗參考和技術(shù)支撐。

參考文獻(References)[1]楊樂平, 朱彥偉, 黃煥. 航天器相對運動軌跡規(guī)劃與控制[M]. 北京:國防工業(yè)出版社, 2010.

YANGLeping,ZHUYanwei,HUANGHuan.Spacecraftrelativemotiontrajectoryplanningandcontrol[M].Bejing:NationalDefenseIndustryPress, 2010. (inChinese)

[2]岳曉奎, 苑云霞. 橢圓軌道相對動力學狀態(tài)轉(zhuǎn)移矩陣[J]. 中國空間科學技術(shù), 2011, 31(1): 42-47.

YUEXiaokui,YUANYunxia.Transfermatrixforrelativedynamicsinellipticorbit[J].ChineseSpaceScienceandTechnology, 2011, 31(1)： 42-47. (inChinese)

[3]SchweighartSA,SedwickRJ.High-fidelitylinearizedJ2modelforsatelliteformationflight[J].JournalofGuidance,Control,andDynamics, 2002, 25(6): 1073-1080.

[4]PluymJP,DamarenCJ.SecondorderrelativemotionmodelforspacecraftunderJ2perturbations[C]//ProceedingsofAIAA/AASAstrodynamicsSpecialistConferenceandExhibit, 2006.

[5]ChenWY,JingWX.DifferentialequationsofrelativemotionundertheinfluenceofJ2perturbationandairdrag[C]//ProceedingsofAIAASpaceConference&Exposition, 2010.

[6]朱仁璋. 航天器交會對接技術(shù)[M]. 北京: 國防工業(yè)出版社, 2007.

ZHURenzhang.Spacecraftrendezvousanddockingtechnology[M].Beijing:NationalDefenseIndustryPress, 2007. (inChinese)

[7]楊嘉墀. 航天器軌道動力學與控制[M]. 北京: 中國宇航出版社, 1995

YANGJiachi.Spacecraftorbitdynamicsandcontrol[M].Beijing:ChinaAstronauticPublishingHouse, 1995. (inChinese)

High-accuracy state propagation of non-cooperative relative orbit in space

ZHU Zhenglong, YAN Ye, YANG Yueneng

(CollegeofAerospaceScienceandEngineering,NationalUniversityofDefenseTechnology,Changsha410073,China)

Abstract：Inordertoresolvetheproblemofrelativenavigationduringthecourseapproachinganon-cooperativespacetarget,ahighlyaccuratestatepropagationmodelofnon-cooperativerelativeorbitmotionwasproposedonthebasisoftherelativeorbitdynamicsequationsandthesecondorderRunge-Kuttaformula.TakingtheJ2perturbationintoaccountandexpandingthetargetorbitequationsintosecondorderTaylorseriesbesidethereferenceorbit,theangularvelocityandtheangularaccelerationofreferenceorbitwerefurtherdeducedandtherelativeorbitdynamicsdifferentialequationwasbuilt.Aftergivingtheinitialvalue,thesecondorderRunge-Kuttamethodwasusedtopropagatetherelativeorbitmotion.Withtheuseofnumericalintegralmethod,themodelcanbewidelyusedwithoutanylimitedconditioninrelativeorbitdynamicsmodelforms;withtheuseofsecondorderRunge-Kuttamethod,thecomputingefficiencyisimprovedwhilethecalculationaccuracyisguaranteed.Twosimulationscenarios,alowEarthorbitscenarioandahighEarthorbitscenario,weredesignedtotestifythegeneralityandprecision.

Keywords：relativeorbitdynamics;J2perturbation;relativeorbitpropagation

doi:10.11887/j.cn.201603014

收稿日期：2015-03-07

基金項目：國家自然科學基金資助項目(11502288)

作者簡介：朱正龍(1987—)，男，湖北宜昌人，博士研究生，E-mail:long420521@163.com； 閆野(通信作者)，男，教授，博士，博士生導師，E-mail:yynudt@126.com

中圖分類號：V412.4

文獻標志碼：A

文章編號：1001-2486(2016)03-081-07

http://journal.nudt.edu.cn