高精度空间非合作式相对轨道状态预报.DOC

1、1高精度空间非合作式相对轨道状态预报 朱正龙，闫野，杨跃能（国防科技大学 航天科学与工程学院，湖南 长沙 410073）摘 要：本文提出了一种基于相对轨道动力学方程和二阶龙格库塔积分公式的高精度空间非合作目标相对轨道预报模型。考虑地球 J2 摄动加速度，将目标轨道方程在参考轨道附近展开，保留至引力加速度差的二阶展开项，并进一步推导考虑 J2 摄动的参考轨道角速度和角加速度，建立相对轨道动力学微分方程；在给定相对轨道初值的情况下，采用二阶龙格库塔积分公式进行相对轨道预报。该模型采用数值积分方法，对相对轨道动力学模型形式没有限制，通用性好，适用范围广；选择低阶龙格库塔公式积分预报，既减小了计算量又

2、保证了计算精度。设置了高轨和低轨两种相对运动仿真场景，仿真结果表明了模型的通用性和精确性。关键词：相对轨道动力学、J2 摄动、相对轨道预报中图分类号：V412.4 文献标志码：A 文章编号： High Accurate State Propagation of Non-Cooperate Relative Orbit in SpaceZHU Zhenglong, YAN Ye, Yang Yueneng(College of Aerospace Science and Engineering, National University of Defense Technology, Changsh

3、a 410073, China)Abstract：Based on relative orbit dynamics equations and second order Runge-Kutta method, a highly accurate state propagation model of non-cooperative relative orbit motion is proposed. Taking the J2 perturbation into account and expanding the target orbit equations into second order

4、Taylor series beside the reference orbit, the relative orbit dynamics equations is built. Then the second order Runge-Kutta method is used to propagate the relative orbit motion. With the use of second order Runge-Kutta method, the computing efficiency is improved while calculation accuracy is guara

5、nteed. Whats more, the model can be widely used without any limiting condition because of numerical integral method is used. Two simulation scenarios, a low Earth orbit scenario and a high Earth orbit scenario, are designed to testify the generality and precision.Keywords：relative orbit dynamics, J2

6、 perturbation, relative orbit propagation收稿日期：2015-XX-XX基金项目：国家 863 计划资助项目作者简介：朱正龙（1987-），湖北宜昌人，博士研究生，E-mail:闫野（通信作者），男，博士，博士生导师，E-mail:空间相对轨道状态预报技术是完成空间接近操作任务的关键技术之一。针对空间非合作目标的相对轨道预报难度更大，不仅需要追踪航天器对非合作目标进行跟踪测量，还需要不依赖目标轨道参数的相对轨道预报模型提供理论支撑。经典相对轨道预报模型，如 CW 模型和TH 模型 12广泛应用于空间交会对接、编队飞行以及其它在轨操作接近任务中，其优点是可

7、获得解析解，使用方便，但是该模型忽略了地球非球形引力摄动，并对相对运动方程进行了线性化处理，适用范围有限和预报精度受到限制，难以满足高精度空间操作任务的需求。Schweighart 和 Sedwick3考虑 J2 摄动的影响，并引入平均参考轨道的概念，建立了线性化的相对轨道模型，适用于近圆轨道的空间编队任务。该模型虽然也可以获得解析解，但是参考轨道的平均化、相对运动方程的线性化处理以及近圆假设均限制了模型的适用范围。Pluym和 Damaren4建立了考虑到中心引力 3 阶、J2摄动引力 2 阶的相对运动模型，该模型仅给出微分方程，并且采用了平均轨道角速度。Chen和 Jing5基于拉格朗日动

8、力学理论，建立了考虑 J2 摄动和大气阻力摄动的精确相对运动微分方程，并给出了相应的精确轨道角速度矢量和轨道角加速度矢量。现有的相对运动模型存在的问题是适用范围受限、通用性差或者难以求解，这都不利于在实际工程应用。本文以接近的空间非合作目标为背景，研究远距离、高速相对运动条件下的高精度相对轨道预报模型。考虑地球 J2 摄动，建立完整相对轨道动力学方程，然后将相对运动方程在参考轨道附近展开，保留主要引力加速度项，并借鉴数值积分的思想，建立模型精度高、计算速度快、通用性好的相对轨道预报解析模型。21 相对轨道动力学1.1 一般相对动力学方程为描述问题的方便，首先定义轨道坐标系：原点 O 位于航天器

9、的质心， ox 轴由地心指向航天器的质心，oy 在轨道面内垂直于 ox 轴并指向运动方向， ox、oy 和 oz 构成右手直角坐标系。设追踪航天器和目标航天器的位置矢量分别为 和 ，追踪航天器的轨道角速度矢量为cRt。定义两个航天器的相对位置矢量为o，根据动坐标系下的矢量微分法则，tc在追踪航天器的轨道坐标系中：* 2ooodt+MERGEFORMAT (1)式中， 表示绝对导数，而 和 表示动坐2dt 标系下的当地导数。假设航天器未受控制力作用，根据牛顿运动定律：* 2cpdtgMERGEFORMAT (2)式中， 表示中心引力加速度差， 为摄动cgpg加速度差之和。联立式* MERGEFO

10、RMAT (1)和式* MERGEFORMAT (2)可得追踪器轨道坐标系下的目标相对运动方程：* 2ooocpgMERGEFORMAT (3)需要注意，采用不同的参考轨道模型，轨道角速度矢量 和轨道角加速度矢量 的表oo达式不同。1.2 精确二体相对运动模型在式* MERGEFORMAT (3)中，忽略摄动加速度和控制加速度，可得到精确的二体相对运动方程：* 2ooocgMERGEFORMAT (4)将 和 投影到追踪器轨道坐标系：,* ,TTTxyzxyzxyz=MERGEFORMAT (5)在二体运动条件下，追踪器轨道角速度矢量和轨道角加速度矢量可表示为 1* 0,0,TTooMERGE

11、FORMAT (6)其中， 为真近点角。目标器和追踪器的中心引力加速度在追踪器轨道坐标系中的投影分别为* ,3,TctctrxyzgMERGEFORMAT (9)* ,21,0TccrMERGEFORMAT (10)中心引力加速度差的表达式为* 320cctcrxyzgMERGEFORMAT (11)其中， 。将式* 22tcrxyMERGEFORMAT (5)至式* MERGEFORMAT (11)带入式* MERGEFORMAT (4)，并展开成坐标分量形式：* 22333/cctttxyrxrxyzrMERGEFORMAT (12)式* MERGEFORMAT (12)即为二体条件下目标

12、相对追踪器的精确运动方程，为时变系数的非线性微分方程组。1.3 精确 J2 摄动相对运动模型文献5在推导相对轨道微分方程时，考虑了 J2 摄动和大气阻力摄动，本文仅考虑 J2 项地球扁率摄动，忽略其它摄动加速度和控制加速度，建立相对运动模型：* 2ooocJgMERGEFORMAT (13)此时，目标器与追踪器的轨道均受到 J2 摄动力的影响。根据文献6，摄动加速度在轨道坐标系中可表示为：* 222413sin, iEJ uiJRriurigMERGEFORMAT (14)在追踪器轨道坐标系中，地球扁率摄动加速度差的计算表达式为 222,TctJioiJttJcriuriugMgg* MERG

13、EFORMAT (15)式中， 和 分别为追踪器和目标器的轨cioti道坐标系相对惯性坐标系的坐标转换矩阵，与瞬时轨道根数有关。轨道角速度矢量 的计算过,Toxyz3程及表达式为* 33200oiuiMMERGEFORMAT (16)式中， 表示绕坐标轴的初等转1,2j换矩阵。进一步整理上式可得：* cosinsicouuiMERGEFORMAT (17)轨道根数的导数可根据高斯摄动方程获得参见文献7，将式* MERGEFORMAT (14)代入高斯摄动方程可得 5：* 23sincoJuiAhrMERGEFORMAT (18)* 23iJiMERGEFORMAT (19)* 223cosin

14、JAuhurrMERGEFORMAT (20)* 23iJhMERGEFORMAT (21)其中， ，21pae为常数。进一步，追踪器轨道223/JEAR角速度为* 323/22sin/0/1cosxJccyzccuihrhreMERGEFORMAT (22)对式* MERGEFORMAT (22)求微分，可求得轨道角加速度 ：2 22465253sinsinsicos0sinxJcJcccyJccz cAuuurhhrrArr * MERGEFORMAT (23)其中， 为追踪器地心距的变化率，且满足c。crv将式* MERGEFORMAT (11)、* MERGEFORMAT (15)、*

15、 MERGEFORMAT (22)和式* MERGEFORMAT (23)代入式* MERGEFORMAT (13)即得到目标器相对追踪器的轨道运动模型。2 相对运动模型的展开从式* MERGEFORMAT (15)可以看出，在计算 J2 摄动加速度差时，还需要获取目标航天器的轨道参数，而对非合作空间目标，其轨道参数常常无法获取。为此考虑将目标航天器轨道在参考轨道航天器附近进行泰勒级数展开。2.1 中心引力加速度的展开根据式* MERGEFORMAT (9)和式* MERGEFORMAT (10)，将 在 附近级数,otg,c展开可得：* ,3, 1102!ototcTiotTi HOgg+e

16、MERGEFORMAT (24)式中， 为三维向量空间中的标准基向量，i表示目标引力加速度 关于 的高阶1HOT,otg展开项。在式* MERGEFORMAT (24)中，作如下符号定义：* ,10otcrPAMERGEFORMAT (25)* ,1 ,12,3TioticrigQMERGEFORMAT (26)经计算推导，可得：* 13201crPMERGEFORMAT (27)* 14201crQMERGEFORMAT (28)* 124031crMERGEFORMAT (29)4* 134012crQMERGEFORMAT (30)根据式* MERGEFORMAT (24)至* MERG

17、EFORMAT (30)，经推导计算，中心引力加速度差的一阶项和二阶项表达式分别为* 13, 2ostc TcxryzgPMERGEFORMAT (31)223,2141Tondiici crxyzeQ* MERGEFORMAT (32)2.2 J2 摄动加速度的展开根据式* MERGEFORMAT (14)，将目标J2 摄动加速度 在追踪器轨道坐标系中的计2,Jtg算式为* 2 2, ,TctJioiJttt riuMMERGEFORMAT (33)为方便计算，将 转化为显含 的函数表达2,Jtg式。根据文献1，目标器在惯性坐标系下的 J2项摄动加速度可表示为* 22225213I et Z

18、XRJYZJRMERGEFORMAT (34)将 投影到追踪器轨道坐标系的关,TtXYZR系式为* cToirxyZzMMERGEFORMAT (35)坐标转换矩阵 的表达式为coi* 313iMERGEFORMAT (36)利用 ，将式* 22,TIttioJttriuJRgMERGEFORMAT (33)化为* 2,2cIiJottMRMERGEFORMAT (37)将式* MERGEFORMAT (34)至式* MERGEFORMAT (36)带入上式，可将 表示2,Jtg为显含相对位置矢量 的函数 ，将2,Jt在 附近级数展开可得：2,Jtg2,Jc* 2222,3, 210Jtt T

19、iJtTi HOg+eMERGEFORMAT (38)式中， 为三维向量空间中的标准基向量，i表示目标引力加速度 关于 的高2HOT2,Jtg阶展开项。作如下符号定义：* 2,20,JtcriuPAMERGEFORMAT (39)* 2,2, ,1,3TiJticrigQMERGEFORMAT (40)记 ， 为任意的角变sin,osc量，经过推导计算可得 4：* 252222231,44735Ecui ui uiii iuiuiic uJRrsssccri PMERGEFORMAT (41)* 2216222254475Ecui ui uiii iiiJRrsssc QMERGEFORMAT

20、 (42)* 226222222547473Ecuiiuiiuii i iJRrscscs QMERGEFORMAT (43)* 22362222254754753EcuiiuiuiiiiiJRrsscs QMERGEFORMAT (44)5根据式* MERGEFORMAT (38)，保留至 J2 摄动加速度差的二阶近似表达式为* 23,112, ,JcicTiruiruQgP+eMERGEFORMAT (45)将上式中的一阶项和二阶项分别记为 和2,1Jstg，即2,Jndg* 2,1,JstcruigPMERGEFORMAT (46)* 23,12,JndicTiQeMERGEFORMAT

21、 (47)若取中心引力和 J2 摄动加速度差的一阶展开项，根据式* MERGEFORMAT (13)可得线性化的相对运动方程：* 2,1,1ostJstgMERGEFORMAT (48)其中， 和 分别由式* MERGEFORMAT (22)和式* MERGEFORMAT (23)给出， 和,1ostg分别由式* MERGEFORMAT (31)和式* 2,1JstgMERGEFORMAT (46)给出。取状态变量为 ，式* ,Tx=MERGEFORMAT (48)写成状态空间的形式：* ttAMERGEFORMAT (49)其中，时变矩阵 的表达式为t* 33212t OIAPMERGEFO

22、RMAT (50)方程* MERGEFORMAT (13)给出了考虑J2 地球扁率摄动的精确相对运动方程，其中 J2摄动加速度差的计算包含了目标航天器的轨道根数，这对于非合作目标的相对轨道计算很不方便。为此，考虑取全部中心引力加速度差和J2 摄动加速度差的前两阶展开项，得到高精度的非线性相对运动方程为* 2,1oJgMERGEFORMAT (51)式中， 和 分别由式* MERGEFORMAT (22)和式* MERGEFORMAT (23)给出， 和og分别由式* MERGEFORMAT (11)和式* 2,1JgMERGEFORMAT (45)给出。3 相对轨道预报模型本节在相对轨道运动微

23、分方程* MERGEFORMAT (51)的基础上，研究快速、高精度的空间非合作目标相对轨道状态预报方法。在进行相对轨道预报时，一般根据当前时刻的状态 及追踪器轨道参数 ，计算出经kxc时间之后的状态 ，可表示为如下ht1kx非线性函数的形式：* 1,kktFMERGEFORMAT (52)为实现相对状态的时间转移，现在的任务就是构造一个计算速度快，计算精度损失小的非线性函数 。借助数值积分的方法，采1:kx用低阶的 Runge-Kutta 积分方法实现相对状态的转移。设微分方程形式的相对轨道模型可表示为* ,;ctxfMERGEFORMAT (53)式中， 表示追踪器的轨道根数，在此假设轨c

24、道根数 可由追踪器上的导航设备精确提供，,k并且在单个状态转移周期 内保持不变。t3.1 线性化的相对轨道预报模型线性系统理论已经十分成熟，而且在使用中十分方便。若线性化的相对轨道运动方程可写为* ttxAMERGEFORMAT (54)根据常微分方程的知识，求解上述连续时变系统，可得* 0,tttxxMERGEFORMAT (55)式中， 状态转移矩阵，其计算公式为0,t* 0eptAMERGEFORMAT (56)对于矩阵的指数运算，其定义为 211exp!kttttAI* MERGEFORMAT (57)从而，当状态转移时间 较小时，可得0t到状态转移矩阵的表达式：* 201,tttIA

25、MERGEFORMAT (58)考虑地球中心引力加速度差和 J2 摄动加速度差的一阶展开模型由式* MERGEFORMAT 6(49)给出，相应的系统矩阵 由式* tAMERGEFORMAT (50)计算，阶数根据需要选择。3.2 非线性相对轨道预报模型对于高精度的状态预报问题，状态方程线性化处理损失了模型的精度，为此须对非线性相对轨道运动方程进行求解。相对轨道运动方程方程由式* MERGEFORMAT (51)给出，记，可将式* MERGEFORMAT (51),Txv写成式* MERGEFORMAT (53)的形式，并且：* 2,1,;coJt fvgMERGEFORMAT (59)基于二

26、阶 Runge-Kutta 公式的相对状态转移公式为* 12112,;/kckckthfxMERGEFORMAT (61)二阶 Runge-Kutta 公式为二阶 精度的2Oh积分公式，具体的状态转移精度由状态转移周期 决定。t4 仿真分析本节从积分方法、积分步长以及相对动力学方程三个方面分析本文建立的相对轨道预报模型的精度。在仿真分析时，采用 STK-HPOP模块生成目标航天器和追踪航天器的轨道数据，并转换为相对轨道数据，作为真实轨道数据用于对比分析模型精度。在高轨和低轨分别设计一个目标器与追踪器先接近后远离的相对运动场景，最小相对距离为 20km，仿真时间为以目标器、追踪器最近距离为中点的

27、前后共 600s 的时间。低轨仿真场景：目标器与追踪器轨道分别运行于 800km 高度的太阳同步轨道上，两轨道面的夹角为 7.69目标器和追踪器先后过轨道交点，最小相对距离约为 20km。高轨仿真场景：目标器运行于地球同步轨道上，追踪器运行于周期为 12 小时、远地点高度低于目标轨道 20km 的椭圆轨道上 。4.1 积分方法对预报精度的影响采用基于二阶 Runge-Kutta 积分公式的相对轨道预报模型，分别在低轨和高轨两个不同场景条件下，仿真分析考虑不同展开项引力加速度相对运动模型的预报精度。图 1 给出了分别基于 Euler 方法和二阶至四阶 Runge-Kutta 方法进行轨道预报的误

28、差曲线，仿真计算步长为 0.1s。仿真结果表明：采用二阶 Runge-Kutta 方法既保证了计算精度，也提高了计算效率。（a ）低轨算例(a) Low Earth orbit example （b）高轨算例(b) High Earth orbit example图 1 采用不同积分方法时的相对轨道预报误差曲线Fig. 1 Relative orbit propagation error with different integral methods4.2 积分步长对预报精度的影响将积分步长分别取为 0.1s、1s、5s 和 10s四种情况分别进行轨道预报，经过 600 预报后，预报误差如表

29、1 和表 2 所示。表 1 和表 2 的结果表明：若采用二阶以上（含二阶）积分方法，模型中计算步长在 0.1s-10s 之间时的误差发散情况相当；而采用 Euler方法时，误差近似呈线性发散，计算精度很低。表 1 采用不同积分步长时的低轨相对轨道预报误差Tab. 1 Relative orbit propagation errors with different integral step in low Earth orbit积分步长/s 0.1 1 5 10Euler 11.08 107.5 532.9 1057RK2 1.2438 1.3182 3.6207 11.572RK3 1.243

30、1 1.2515 1.2962 1.3683RK4 1.2431 1.2515 1.2966 1.3720表 2 采用不同积分步长时的高轨相对轨道预报误差Tab. 2 Relative orbit propagation errors with different 7integral step in high Earth orbit积分步长/s 0.1 1 5 10Euler 3.113 30.89 154.4 308.8RK2 0.3687 0.3688 0.3730 0.4031RK3 0.3678 0.3688 0.3710 0.3784RK4 0.3678 0.3688 0.3710

31、0.37844.3 相对动力学方程对预报精度的影响当相对轨道动力学方程保留不同阶数的地球引力加速度差时，得到的相对轨道预报模型精度不同。图 2（a）和图 2（b）分别为低轨和高轨场景下的相对轨道预报误差曲线。仿真结果表明：在 600 秒的相对运动过程中，二体相对模型TB及其一阶展开式TB（1）、二阶展开式 TB（1,2）在低轨预报误差均在 50m 以上，而在高轨的预报误差小于 2m；考虑 J2 摄动及其一阶、二阶展开项的模型，即TB+J2、TB+J2(1)和TB+J2(1,2)，在低轨的预报误差均在 10m以内，而高轨预报误差在 1m 以内；此外两个线性模型，即TB(1)和TB(1)+J2(1

32、)在高轨的预报误差在 2m 以内，在低轨的预报误差均超过 1000m。综上可知：在高轨，采用二体模型的一阶展开式（即著名的 Lawden 方程）建立预报模型即可获得较高的预报精度，预报误差在 2m以内；在低轨，除考虑中心引力差的一阶展开项，还必须考虑中心引力的二阶展开项以及 J2摄动加速度差的一阶展开项，才能获得预报误差小于 10m 的模型，若进一步考虑 J2 摄动加速度差的二阶展开项，600 秒内的预报误差可缩小到 2m 以内。J2 摄动是最主要的轨道摄动因素之一，考虑到 J2 摄动加速度的大小近似与航天器地心距的 5 次方成反比，而高轨航天器的地心距约为低轨的 6 倍，由此可知高轨 J2

33、摄动加速度的影响很小，而低轨影响明显。此外，由于相对距离远，摄动加速度差的高阶展开项不再是小量，因而忽略会带来较大的预报误差。（a ）低轨算例(a) Low Earth orbit example（b）高轨算例(b) High Earth orbit example图 2 采用不同动力学方程时的相对轨道预报误差曲线Fig. 2 Relative orbit propagation errors with different equations of relative orbit dynamics4 结论空间非合作目标相对轨道预报是实现空间操作的一项重要技术环节。本文考虑主要轨道摄动因素，并基于

34、二阶 Runge-Kutta 方法，综合考虑计算精度和计算效率，提出的模型原理简单，通用性好，方法简单，计算速度快，预报精度高，可以实现对高速相对运动空间非合作目标的相对轨道预报，在给定的时间内，其预报误差在米量级。本文提出的方法对解决未来更高精度的空间操作任务能够提供经验参考和技术支撑。参考文献1 杨乐平, 朱彦伟, 黄焕. 航天器相对运动轨迹规划与控制M. 北京:国防工业出版社 , 2010.Yang Leping, Zhu Yanwei, Huang Huan. Sapcecraft relative motion trajectory planning and control M. B

35、ejing: National Defense Industry Press, 2010. (In Chinese)2 岳晓奎, 苑云霞. 椭圆轨道相对动力学状态转移矩阵J. 中国空间科学技术, 2011, 31(1): 42-47.Yue Xiaokui, Yuan Yunxia. Transfer Matrix for Relative Dynamics in Elliptic Orbit J. Chinese Space Science and Technology, 2011, 31(1), 42-47. (In Chinese)3 Schweighart S A, Sedwick R

36、 J. High-Fidelity Linearized J2Model for Satellite Formation Flight J. Journal of Guidance, Control, and Dynamics, 2002, 25(6): 1073-1080.4 Pluym J P, Damaren C J. Second order relative motion model for spacecraft under J2 perturbations C. AIAA/AAS Astrodynamics Specialist 7Conference and Exhibit, K

37、eystone, Colorado, 2006.5 Chen Weiyu, Jing Wuxing. Differential equations of relative motion under the influence of J2 perturbation and air drag C. AIAA Space 2010 Conference & Exposition, Anaheim, California, 2010.6 朱仁璋. 航天器交会对接技术M. 北京: 国防工业出版社, 2007.Zhu Renzhang. Spacecraft Rendezvous and Docking Technology M. Beijing: National Defense Industry Press, 2007. (In Chinese) 7 杨嘉墀. 航天器轨道动力学与控制M. 北京: 中国宇航出版社, 1995Yang Jiachi. Spacecraft Orbit Dynamics and Control M. Beijing: China Astronautic Publishing House, 1995. (In Chinese)