1地月平动点拟周期轨道模型及其特征
1．1坐标系定义首先定义文中使用的3个坐标系，其示意图如图1所示。J2OOO地心赤道惯性坐标系(J2ooogeocentricequatorialcoordinate，J2000)0一XY2：其定义参见文献E23。地心旋转坐标系(geocentricrotatingcoordinate，GRC)O—xy2：原点在地球质心，2轴在地月连线方向指向月心，轴指向瞬时月球绕地球运动的轨道角动量方向，轴与其他两个轴构成右手直角坐标系。L2点瞬时旋转坐标系(L2rotatingcoordinate，LRC)I2一2222：原点在地月I2点，22、2、Y2轴分别与2、、轴平行。3个坐标系之间的转换关系可参见文献E23。
1．2高精度动力学模型以及圆形限制性三体模型高精度的星历模型对于计算精确地月平动点拟周期轨道是必不可少的。文中使用了建立在GRC中的高精度星历模型。文献[15]已经对该模型的精度进行了验证并证明了该模型在数学上完全等价于J2ooo惯性系下的星历模型。文献[161则将该模型求得的轨道在STK软件中进行了验证，结果表明该模型精度与STK软件中模型精度相当。现将模型的推导过程简单描述如下。如若忽略其他行星的引力影响以及地球的非圆形摄动影响，根据牛顿经典力学原理，地月平动点探测器在J2OOO坐标系中的动力学模型可表示为一一U!E，4-c一rM+,rp一一”—一十Us(—一一(1)式中，rP为探测器到地球质心的位置矢量；rM，rs分别为月球和太阳相对地球质心的位置矢量；T'pM为探测器与月球之间的距离，可以表示为lr一rPll；rPS为探测器与太阳之间的距离，可以表示为IIR一r尸l。由于J2OOO为惯性坐标系，GRC在空间运动的绝对角速度∞即为其相对于J2OOO的角速度。记探测器在GRC中的相对速度为、相对加速度为，则得到一群-42~o×+∞×(∞×r尸)-4西×r尸(2)联立式(2)和式(1)，得到一一2∞×一∞×(∞×rP)一面×rP一rP+r；一一4-Ps。c一一舞㈤对式(3)进行简单的分析，含有的项是太阳对于探测器的直接弓1力，这一项直接作用在探测器上。而太阳的间接引力以及月球偏心率的影响则体现在国的表达式中。需要指出的是，由于GRC的角速度即为月球绕地球的旋转角速度，因而∞的表达式由月球的运动方程决定。参考文献[15]以及一般摄动方法，不难推导出角加速度在GRC中的表示方式如下：通过悄个推导过程可以看出，式(s)中涉及的有关星历的运算都足矢狱的点乘，与坐标系无关，因而可以使用X2000.}，标准?}}.历数据，例如DE4210进一步分析，如果忽略式(5)中太阳直接引力对探测器的影响、太阳间接引力对月球轨道的影响以及月球轨道本身的偏心率，即认为(fps-0>},=0,y=0,ree=0)，则式(i)简化为圆形限制性三体模型由推导的过程可知，式(5)理论上与绝对动力学模型完全等价，而圆形限制三体问题模型则是文中高精度星历模型的一种特例情况。
1.3拟周期轨道特征以STK}ASYYO},TIYOY中的地月平动点拟周期Halo轨道初值所生成的轨道为例，对拟周期轨道特征加以说明。该轨道起始于2018年12月30日]2:(>0(格林尼治时间)，持续时问为39.29天(约为2个半周期)，最大直径约为10()。()kmo图2(a)是该轨道在GRl'坐标系中的图像，显然，从地心上来看，轨道已经不再呈现周期性。而图2（b)则是在瞬时的LRC上看到的同一条拟周期Hal。轨道。在I.:点上看，同样的轨道呈现出一定的周期性。结果与文献巨77一18}的结论是一样的。与日地系统相比，地月系统平动点轨道的周期更短、第二主天体的轨道偏心率更大、太阳作为第三引力体对地月系统的引力影响更强烈，圆形限制性三体问题中的周期轨道并不存在于地月平动点附近，取而代之的是一种“拟”周期轨道的形态。但是在瞬时平动点旋转坐标系中，该轨道依然作常地接近于周期轨道。
2平动点拟周期轨道设计方法
根据引言中对拟周期轨道设计问题的分析可知，Halo轨道的周期性在高精度星历模型中并不存在，无法使用传统的微分修i1:来设计轨道。除此之外，圆形限制性三休问题的解析解y真实值相差较大，无法为高精度星历模型提供理思的迭代初仇。因此，想要通过多步打靶法求解高精度星历模}g条件下的轨道则必须找出适用于高精度星历模J}'J的新拼接点信息。针对一上述问题，本文提出了一种基于地月平动点拟周期轨道特征的轨道设计方法，其具体设计思路如图:3所示。联系第I.a-}i中所得结沦，由于真实的地球一月球一飞行器的系统会受到月球偏心率以及太阳直接引力和间接引力的影响在这种条件下，原本圆形限制性三体条件下规律的轨道形态会发生变化。由于月球偏心率以及太阳间接弓}力的影响，I点与地心之间的距离并不是常值(即I.:点的位置不是固定的)，因此，在地心旋转坐标系下，轨道会呈现沿着地月连线“见动”的情况。虽然.在地心旋转系下，其形态的规律性已经不存在，但是在瞬时平动点坐标系下，轨道的形态与圆形限制性二体条件下轨道的形态仍然非常相似通过I几述分析可以找出地月平动点拟周期轨道在高精度星历模型条件下，其形态特征中的不变流形。除去不变的部分，其形态特征中变化的部分是由于I:点在地月空间的位置不固定肉造成的。由此，对f拟周期轨道形态特征中不变的部分，可以通过圆形限制性二体条件下平动点旋转坐标系下的轨道来近似而对于其形态特征中变化的部分.则可以通过求解瞬时f点在地心旋转坐标系下的状态来表小。山此川一以认为.在选定轨道设计的历元时间后，在圆形限制性刃体条件下的连续轨道上取一系列已知时问间附的点作为拼接点.按照对应时刻条件下的瞬时I.._点状态数据转换到地心旋转坐标系中。转换后的点将I卜常接近丁地心旋转坐标系中一条满足精确动力学条件的轨道I几的点。这是由于1.,点状态数据包含有月球偏心率以及太阳问接引力的影响。通过对圆形限制性三体条件下已知时问间隔的点叠加I:点状态信息，可以找出含有摄动因素的拼接点。由于拼接点是基于近似轨迹上的点经过坐标转换求得的，转换后的点在精确的力学模型中积分后无法形成一条封闭曲线，因此，需要用多步打靶法来实现轨道的拼接轨道设计的步骤描述如下。步骤1生成圆形限制性三体条件下轨道。使用微分修正生成圆形限制性三体条件下I:点旋转坐标系下的连续Halo轨道。步骤2拼接点求取。在连续Halo轨道r.选取拼接点并选取任务轨道的历元时间T.。标记各个拼接点为N,.其状态为X，对应时刻为T,(i=1.2.""..rn)a求解各个于时刻所对应的瞬时I.:点状态，通过矢量之间几何关系.将瞬时I.:点状态叠加到拼接点的信息中。步骤3轨迹拼接。将拼接点N，在地心精确动力学模型中由T时刻积分到T，十1时刻，形成，。一1段弧段轨道。使用变时间的多步打靶法对轨道进行拼接。各步骤中所含具体内容将在第2.1节一第2.3节中详细说明。
2.1圆形限制性三体条件下连续Halo轨道设计
圆形限制性三体条件下的Hal。轨道设计主要是利用周期轨道本身存在的对称性以及微分修正方法对I2ichardson提出的三阶解析解进行修正，从而得到三体条件下精确轨道初值.再由精确初值积分得到轨道。具体表述为:将式(6)的原点移到姚点，得到LRC系下的圆形限制性三体动力学方程.并结合摄动理论的思想，运用Lindstedt-Poincare法，得到圆形限制性三体条件下的Hal。轨道在LRC'中的三阶近似解.如下式所示:通过改变:’的值，可以得到周期轨道的三阶近似解，式中其他参数的含义及取优参考文献在利用二阶解析解得到初值后，根据Halo轨道关于r一二平面对称的特征，以周期轨道上垂直穿过‘二一二平面的点于作为迭代初值，并以该轨道再次穿越二一二平面的时刻作为终端状态。引为对于初始向量的修正量。为简化修正过程，固定，得到修正关系式中，电为状态转移矩阵中中的元素。通常迭代4^}6次就可以得到精度较高的结果。具体求解过程，参见文献[6}0由此，得到LRC中圆形限制性三体条件下的连续轨道。
2.2考虑摄动条件的拼接点求解
2.2.1拼接点的位置选取
为了方便表示拼接点在Hal。轨道上的位置，首先给出有关定义Halo轨道相角:的定义。:是一种非几何概念下的角度，类似于真近点角的概念，在文中只是用来表示拼接点在Hal。轨道上的位置。相角:的起点00选在轨道垂直穿越二一2平面朝向y轴正方向继续延伸的位置，角度增加方向与轨道延伸方向一致。相角的定义为(9)式中，t表示轨道上的任意时刻;t。表示轨道穿越x-2平面朝向y轴正方向继续延伸的时刻，即00的时刻;p表示Hal。轨道的周期。图4给出投影在二一Y平面内，一个周期Halo轨道相角的位置。由图4可以看出，虽然相角二之间相差均为900，但是实际其位置分布并非是等分的。由文献[12]中的分析可知，Hal。轨道动力学敏感性随轨道曲率增加而增加。因而，需要在轨道曲率较大的位置设置较密集的拼接点，过少的拼接点可能导致仿真的发散。对于文中第2.1节中求解的连续轨道，对其每一周期分别选取相角:为450,1350,2250以及3150的点作为拼接点。需要指出的是.拼接点的位置无需固定为上述4个位置，只要保证在轨道曲率较大时较密集即可。
2.2.2叠加瞬时平动点的状态信息
由于三体模型所用的参数是一些平均值，对于历元时间没有要求，但是，带有星历的模型必须严格依据具体时间(格林尼治时间)，不同时刻对应不同的参数条件，因而在坐标转换之前必须确定历元时间。每一个拼接点的信息必须包含时间和该时间对应的位置速度状态量，例如:N;（T;,X;)。需要指出，此时X，的信息是在LRC中，需要通过矢量关系将瞬时L2点的状态信息添加到拼接点的信息中。图5表示出了矢量之间的几何关系，R为地心到探测器的矢量在GRC中的表示;斌.。为L2点在GRC中的位置矢量，其导数表示L:点在GRC中的速度矢量;r表示L2点到探测器在GRC中的位置矢量。式中，RM和V,',}分别是月球的位置速度在J2000'}性坐标系中的表示，由于向量的点乘与坐标系统无关，因此，可以使用J200。中的星历数据直接计算，而r“则正是已知的拼接点在LRC下的信息戈。
2.2.3参数几的求取
注意，通常情况下，对于又的取值可以选择只与地月的质量比有关而与历元时间无关的经典共线拉格朗日参数[19]。当取地月质量比为81.30065597时，得到几一1.16783268238542。但是，由于月球轨道的偏心率以及太阳引力的影响，共线拉格朗日参数并不是固定的，若要实现精确的轨道求解，可以利用迭代的方法求得轨道拼接弧段时间内久的精确解，具体求解方法如下。对于参数几的求取就是对太阳引力摄动情况下L2点的精确位置的求取。因而，设式(5)以[、(‘。)，y(t9),2(tp),vs(tp),vv(to),}2(tp}}Tr巨nYM(tp),O,O,aYm(to),O,O}T为初始条件，其解具有如下形式:P=[二.y，二，'，夕，'」T}-,y,2}其中从为常数。当久取地月质量比时，点P应在平动点L=以rM,O,叼T附近。根据上述分析，可以求取使得式(13)达到极小值的参数几。J(a,to，△!)一厂}(P一L>CP一L}dr-}(P(}t)一LC}t))‘(P（}t)一LC}tW}}（13}式中从是一个与初始时刻t。和持续时间△t有关的量。然而，对具体的在轨任务而言，t。以及△t是给定的。因此，可以给出J(d,to,0t)取极小值时，久的迭代算法为几i+1一几，一J(},,tp,fit)J(}+}},tp,0t)一J(}l},to,fit)(14)i=0，1，2，…式中，以为小的正数，将△t作为变量，对式(14)求导，得到J(.l,tp,}t)=2(P（t)一L(t;T(P（t)一L(t;+2(P(t)一L(t))丁(户(t)一L(t))t(15)且JUfto,0)=0。联立式(15)与式(5}，对其从t。到to+pt积分得到J(}l.,tofpt)a迭代初值几。可以取经典拉格朗日参数[u9}。给定初始时刻t。以及时间间隔pt。利用式(14)计算得到对应时间区域内的几。
2.3轨迹拼接
将转化到GRC的拼接点按照之前设定好的时间间隔用带有星历的动力学向前积分，积分后的轨道会出现不重合的特征。例如，如图6所示，N,点从T:时刻，以状态戈向前积分到T2时刻，状态成为X}2，表示为NsCT2,Xl2)o而Ta时刻应该对应于拼接点N2点的状态凡。此时在N2点出现了两个不同的状态。同理，以N:点的状态向前积分到T3时刻后成为N3(T3fX23)与N3点也不重合。这时轨道设计问题就变成了典型的轨道拼接问题。斩道拼接根据文献[7]提出的多步打靶法，通过位置拼接和速度拼接两个部分来完成。
2.3.1位置拼接
为了解决轨道拼接问题，我们首先通过修正N〕点速度来调整N2点位置，从而实现与N:点的位置拼接。微分修正关系为式中，巾。是Ni点到N2点的状态转移矩阵。中第i行第J列对应的元素。N，点到N2点的修正完成后，以N2点为起点，N,点为终点，进行第二段修正，直至所有的轨迹在位置上都连续。
2.3.2速度拼接
位置修正完成后，虽然N2点与N2点的位置重合，但是仍然存在速度跳变△}2，因此还要继续进行速度修正，从而使得轨迹在速度上连续。记N从点的位置状态r;一「二，，y，2〕Ty;,2;，速度v;=}x，2}]T，加速度a;a以N,点为初值积分到兀时刻所得到的状态记为r:一，v2一,a2_，并记其状态转移矩阵为梦‘;以N。点为初值逆向积分得到T2时刻所得到的状态记为r2+,V2+,a2十，并记其状态转移矩阵，状态转移矩阵表示为从而可以将拼接点N2处的速度跳变△v2表示成NN2,N33点的位置r和时间T的函数，并通过改变这3点的位置和时间来使得△v2趋向于零。经过推导，可以得到关系式通过式(17)可以得到拼接点从速度的微分修正改变量。式(17)是仅有3个拼接点情况下的修正公式，将此方法类推到，-1段轨迹，二个拼接点，则得到相应的微分修正关系式为式中在完成速度修正后，拼接点的位置仍然会发生改变，所以需要继续对拼接点位置进行修正，如此反复修正，直到位置和速度误差都在允许范围内则停止迭代，即得到连续的拼接轨道。
3数学仿真
仿真以设计一条起始历元时间为2013年10月1日12:00（GMT)，时长为11。天的地月平动点拟周期轨道为例，来具体说明设计的过程。设轨道的最大直径约为10000km，设计步骤如下。步骤1连续Hal。轨道设计根据式((7)中的二阶近似解公式，井令r"=0,n=1,f1:一0.06137(即x方向幅值为5000km),X=1800，可以得到圆形限制性三体轨道初值的近似解。上述参数条件下的Hal。轨道周期为15.7167天。按照式((8)通过微分修正方式得到长度为110天(约7个周期)的连续轨道，如图7所示。由于图7中轨道是利用轨道关于x一二面对称的特点以及微分修正所求的。因而轨道的起始点为轨道穿越x-2平面并朝向y轴正力一向继续延伸的位置，符合相角:=00的定义，由此可知式((9)中t。即为历元时刻。进而可以计算出仿真时间内相角相对于时间的变化，如图8所示。图8中，实线表示仿真时间内相角随时间的变化。点化线分别对应450.1350,225“以及31504个角度。通过仿真可以看出，在仿真时间内轨道28次经过预先设定的拼接点相角，记录这28个点的位置以及时间信息，选为拼接点。从图7中还可以看出，起始点时刻的位置相角为00，终点时刻的位置相角为359.60，不与预先设定的相角重合，因而也需选为拼接点。图7中的轨道共选择拼接点30个，形成29段弧段。接下来，根据第2.2.3节中a的求解方法，分别求解29时间段所对应的几值。按照式(11)，结合拼接点时间信息T所对应的星历数据，将30个拼接点逐一转化到GRC中，并向前积分相对应的时间间隔。仿真形成的不闭合轨迹如图9所示，图中只画出了最初的5段曲线。图9中，实线表示由拼接点积分形成的不闭合弧段。从图9中可以看出，仿真结果与第2节中的分析一致，由于拼接点是基于近似轨迹上的点经过坐标转换完成的，转化后的点在精确的力学模型中积分后无法形成一条封闭曲线，但是非常接近于一条封闭曲线。步骤3轨迹拼接采用第2.3节中多步打靶的方法，将轨道拼接完成，其轨迹如图10(a)所示，图10(b)一图10（d)分别为其在3个平面内的投影。文中所述的方法适用于摄动条件下各类地月平动点拟周期轨道的设计，包括拟周期Lyapunov轨道、拟周期L1S-sa}ous轨道和拟周期Halo轨道。因为摄动条件下，上述轨道在瞬时平动点坐标系中均具有在与圆形限制性三体条件下相类似的形态特征。因此，同样可以使用圆形限制性三体条件下平动点旋转坐标系下的Lyapunov或者Lissajous轨道来近似高精度星历模型条件下Lyapunov或者Lissa-jous轨道特征中不变的部分。通过求解瞬时平动点在地心瞬时旋转坐标系下的状态来表示其形态特征中变化的部分。由于文章篇幅有限，此处不再赘述。
4结论
本文利用地月精确动力学模型研究了设计地月平动点拟周期轨道的方法，得出如下结论。
(1)建立了带有真实星历数据的高精度地月平动点拟周期轨道动力学模型，该模型中考虑了太阳的直接引力、间接引力以及月球偏心率的影响。
(2)利用地心旋转坐标系中的高精度动力学模型研究了采用轨道拼接方法迭代计算拟周期轨道的方法，并详细说明了考虑摄动条件下拼接点的选取方式。
(3)通过数值仿真，分别使用3种类型的拟周期轨道对算法进行了验证，计算了具有30个节点，周期为11。天的拟周期绿色农业论文轨道。仿真结果表明，算法能够有效地设计出所需轨道。
(4)文中所述的方法适用于摄动条件下各类地月平动点拟周期轨道的设计，包括拟周期Lyapunov轨道、拟周期Lissajous轨道和拟周期Halo轨道。
作者：钱寞婧 单位：哈尔滨工业大学航天工程系
