下载:
-
近年来，随着航天技术的快速发展，深空探测越来越受到世界各航天大国与组织的重视。按照我国航天事业发展的指南，未来航天的3个重点方向之一就是深空探测[1]。因此，与深空探测密切相关的各项基础和应用研究便成为了我国航天事业中的一个热点研究方向。
不同于经典二体问题下的近地空间航天活动，深空探测所涉及的运动模型更多的是三体问题。由于动力学本质上的区别，使得三体问题相对于二体问题更加复杂。首先，从数学的角度上，三体问题是一个强非线性时变系统，所以必然存在强非线性系统所具有的共性问题，即状态终值对初值具有极高的敏感性与不确定性，导致三体下的航天器运动具有很强的“混沌”性；其次，在三体问题中，航天器的运动，尤其是在共线平动点附近的运动，极易受到外界环境摄动因素的影响，导致其后的运动呈指数发散趋势，使得运动具有强不稳定性；最后，在三体问题中共有18个自由度，而目前只找到13个运动积分，这使得三体问题下的航天器运动不存在类似二体问题那样的精确解析解[2]。所以，退而求其次，为获得三体问题下运动的近似解析解，通常需要对三体问题进行适当的简化，其中最基本、最常用的模型依次是圆型限制性三体问题（Circular Restricted Three-body Problem，CRTBP）模型和椭圆型限制性三体问题（Elliptic Restricted Three-body Problem，ERTBP）模型。
在CRTBP和ERTBP的研究中，平动点及其附近的周期/拟周期轨道由于其独特的位置优势使其在深空探测任务设计中占据着重要的地位。距离地球最近的三体系统是日–地/月三体系统和地–月三体系统，这两个三体系统的平动点轨道一直是许多已经提出或完成了的深空探测计划的工作轨道，如美国国家航空航天局（National Aeronautics and Space Administration，NASA）的ISEE、NGST计划，欧洲航天局（European Space Agency，ESA）的GAIA、DARWIN计划，以及我国的“夸父”计划、“鹊桥”通信中继卫星计划等。通常情况下，对平动点及其附近运动的研究多依靠数值方法实现，如数值积分、微分修正、数值优化等，然而，数值方法本身具有其固有的劣势，即其本质上是一个试错过程，不能完全体现运动的本质和规律，并且很大程度上依赖于设计者的经验。与数值方法相反，解析方法却能够帮助更加深入地研究平动点附近的运动，探究运动的特性和规律，为后续的深空探测任务轨道设计提供一定的先验信息，因此三体问题下平动点附近运动的解析研究变得日益重要。
CRTBP是最简单的三体运动模型，其没有考虑主天体公转轨道的偏心率以及各种环境摄动力，主要用来研究三体轨道最基本的动力学特征。目前关于CRTBP平动点及其附近运动的解析研究已经较为全面，Gómez[3]研究了CRTBP中平动点位置、平动点附近轨道的一阶近似解析解，如Lissajous轨道、Lyapunov轨道、halo轨道；Richardson[4]利用Lindstedt-Poincaré方法构造了CRTBP下共线平动点halo轨道的三阶近似解析解；Masdemont[5]研究了CRTBP下共线平动点附近运动的双曲形态和中心形态，将不变流形扩展成双曲振幅和中心振幅的幂级数形式；Jorba等[6]同样采用Lindstedt-Poincaré方法并结合Normal Form Scheme方法半解析地构造了CRTBP下共线平动点动力学方程的解；Rithwik等[7]使用基于优化技术的梯度和非梯度方法研究了halo轨道的设计问题；刘刚等[8]以CRTBP中的周期轨道作为迭代初值，用星历表数据对轨道进行拼接获得了所需的拟周期轨道；陶雪峰等[9]在CRTBP框架内，针对不同类型平动点轨道的特点，通过设定轨道特征点，以轨道闭合程度为目标函数，提出了一种基于优化算法的平动点轨道生成方法；Lara等[10]使用摄动方法和归一化方法来代替Lindstedt-Poincaré方法，获得了平动点附近Hill问题的解析解；郑越等[11]建立了一种改进型庞加莱截面图，结合状态转移矩阵和打靶法，提出了一种计算周期轨道的新方法；Lei等[12]构建了CRTBP下三角平动点附近运动的高阶解析解，并详细讨论了解的收敛性；Liang等[13]利用改进的参数变分方法，获得了三角平动点附近运动的高阶近似解析解；Qian等[14]受振动动力学中的非线性模态启发，结合多项式展开和微分修正方法获得了三角平动点附近垂直周期轨道的解析解和数值解；翟冠峤等[15]从模态运动的角度出发，分析了三角平动点附近周期轨道，通过多项式展开法构建出了周期轨道3个运动方向之间的渐近关系。
虽然在一些问题的分析中，CRTBP能够满足应用背景的要求，但是在研究平动点运动稳定性以及航天器在平动点附近运动的控制问题等方面，CRTBP显然无法满足要求[16]。相比之下，ERTBP由于考虑了主天体公转轨道的偏心率，因此能够更准确地描述三体系统内的航天器运动情况，同时也使得对ERTBP的基础研究更具有研究意义和工程应用价值。然而，相比于CRTBP，ERTBP下平动点附近的运动研究更为复杂，主要体现在：①CRTBP对应的是自治系统，而ERTBP对应的是非自治系统；②CRTBP下的动力学方程为定常系统，而ERTBP下的动力学方程为时变系统，运动积分（Jacobi积分）不再存在[17]；③CRTBP中平动点附近存在周期轨道，而在ERTBP中平动点附近不存在周期轨道，只存在相应的拟周期轨道[18]。上述因素直接导致目前有关ERTBP下共线平动点附近运动的研究相对较少，国际上普遍认为对ERTBP的研究还不够深入[19]。
目前有关ERTBP下平动点附近运动的解析研究主要有：雷汉伦[20]推导了ERTBP下的无量纲运动方程，并指出在ERTBP中同样存在5个平动点，也同样分为3个共线平动点和2个三角平动点；Hou等[21]半解析地研究了ERTBP下共线平动点附近的Lissajous轨道和halo轨道的高阶解，并且将该级数解应用于地–月三体系统和日–地/月三体系统；Lei等[22]利用Lindstedt-Poincaré方法，以CRTBP下的一阶近似解析解为初值，构造了脉动会合坐标系下ERTBP共线平动点halo轨道和Lissajous轨道对应的不变流形的高阶解析解；Hou等[23]研究了ERTBP下三角平动点附近的拟周期运动；Lei等[24]同样利用Lindstedt-Poincaré方法构造了ERTBP下三角平动点附近运动的高阶解析解，并将其应用于地球至地月三角平动点周期轨道的低能转移轨道设计中；Gong等[25]在脉动会合坐标系下设计了ERTBP下的太阳帆航天器的周期轨道，并对其稳定性进行了研究，最后利用时变线性二次型调节器实现了周期轨道的稳定控制；Koh等[26]研究了平面ERTBP下的姿轨耦合情况下的运动周期解；Gómez等[27]提出了多点打靶方法，用于求解限制性三体问题下的拟周期轨道。
除此之外，相关文献也对摄动情况下的CRTBP和ERTBP平动点及其附近的运动解析解进行了研究。Narayan等[28]利用基于特征指数的摄动技术研究了考虑主天体辐射压下ERTBP三角平动点附近运动的稳定性；Farquhar等[29]构建了考虑太阳引力和月球轨道的偏心率时的地–月三体系统L2平动点附近的拟周期运动的三阶解析解；Singh等[30]利用Lindstedt-Poincaré方法研究了主天体扁率和辐射压下的平面限制性三体问题共线平动点周围Lyapunov轨道的三阶近似解析解；Hou等[31]解析地研究了JPL星历下真实地月系统共线平动点和三角平动点的拟周期运动；Ansari[32]研究了考虑主天体扁率和辐射压下的圆型限制性四体问题平动点附近的周期轨道；Papadakis[33]通过平面轨道的分岔方法研究了圆型限制性四体问题中的三维对称周期轨道；Jagadish等[34]研究了考虑主天体辐射的摄动限制性三体问题共线平动点附近周期运动的半解析解和数值解；Howell等[35]提出了两层微分修正方法，用来求解真实力学模型下的拟周期轨道。
虽然上述参考文献已经获得了ERTBP下平动点附近周期轨道及对应流形的高阶近似解析解，但均是以CRTBP下的一阶近似解析解（而非ERTBP下的一阶近似解析解）为初值，经过Lindstedt-Poincaré 方法迭代或者打靶法、微分修正等数值方法迭代获得的。然而，如果能够获得ERTBP下平动点附近运动的一阶近似解析解，不仅相比CRTBP下的解析解更具一般性[20]，也能取代CRTBP下的近似解析解，为ERTBP下共线平动点附近运动的高阶解析解研究（如文献[22]）提供一个更加精确的初值，使得迭代次数减少、收敛速度加快，同时也必将有利于进一步完善ERTBP下平动点附近运动的解析研究，为ERTBP的研究奠定更加坚实的基础。基于此，本文参考借鉴CRTBP中的研究方法，对ERTBP共线平动点附近运动的近似解析解展开了相关研究。
-
限制性三体问题研究的是一个小天体或航天器在两个大天体（主天体）引力场中的运动规律，其中小天体或航天器的质量相比于两主天体可忽略不计，不会对两主天体之间的相互运动产生影响。当主天体围绕公共质心作圆轨道运动时，称为CRTBP；当主天体围绕公共质心作椭圆轨道运动时，称为ERTBP，下面分别进行介绍。
-
CRTBP是最简单的三体运动模型，其没有考虑主天体公转轨道的偏心率以及各种环境摄动力，主要用来研究三体轨道最基本的动力学特征。为描述航天器运动的方便和计算的简化，在CRTBP中通常需要对相关物理量进行无量纲化处理，并以时间作为独立变量。相应的质量
$[M]$ 、长度$[L]$ 和时间$[T]$ 的归一化单位取为$$\left\{ \begin{split} & \left[ M \right] = {m_1} + {m_2}\\ & \left[ L \right] = {L_{12}}\\ & \left[ T \right] = {(\dfrac{{{L_{12}}^3}}{{G({m_1} + {m_2})}})^{1/2}} \end{split} \right.$$ (1) 其中：
${m_1}$ 、${m_2}$ 为两主天体质量；${L_{12}}$ 为两主天体之间的距离；$G$ 为万有引力常数。为了更直观、更清晰地描述航天器在CRTBP中的运动，通常选择在会合坐标系（也称为质心旋转坐标系）下进行研究。如图1所示，会合坐标系的原点位于两主天体的公共质心，x轴由较大主天体指向较小主天体，z轴与主天体系统角动量方向平行，y轴满足右手定则[19]。
图 1质心旋转坐标系
$O - XYZ$ 和平动点旋转坐标系${L_2} - xyz$ Figure 1.Barycentric rotating coordinate system
$O - XYZ$ and libration-point centered rotating coordinate system${L_2} - xyz$ 假设
${{x}} = {[x,y,z]^{\rm T}}$ 为航天器在会合坐标系中的位置矢量，则CRTBP下航天器运动的动力学模型为$$\left\{ \begin{split} & \ddot x - 2\dot y = \dfrac{{\partial \varOmega }}{{\partial x}} \\ & \ddot y + 2\dot x = \dfrac{{\partial \varOmega }}{{\partial y}} \\ & \ddot z = \dfrac{{\partial \varOmega }}{{\partial z}} \\ \end{split} \right.$$ (2) 根据文献[36],CRTBP下共线平动点附近运动
$\delta {{x}} = {[\delta x,\delta y,\delta z]^{\rm T}}$ 的近似解析解为$$\left\{ \begin{split} \delta x(t) = & {A_1}{{\rm e}^{{\lambda _1}t}} + {A_2}{\rm e}^{ - {\lambda _1}t} + {A_3}\cos ({\rm Im} ({\lambda _3}))t + \\ & {A_4}\sin ({\rm Im} ({\lambda _3})t) \\ \delta y(t) = & {k_1}{A_1}{{\rm e}^{{\lambda _1}}}t - {k_1}{A_2}{\rm e}^{ - {\lambda _1}t} - {k_2}{A_3}\sin ({\rm Im} ({\lambda _3}))t + \\ & {k_2}{A_4}\cos ({\rm Im} ({\lambda _3})) t\\ \delta z(t) = & {A_5}\cos ({\rm Im} ({\lambda _5})t) + {A_6}\sin ({\rm Im} ({\lambda _5})) t \end{split} \right.$$ (3) 式中，
${A_i}(i = 1, \cdots , 6)$ 是由初始条件确定的积分常数，参数${k_1},{k_2},{\lambda _1},{\lambda _3},{\lambda _5}$ 的详细计算公式可以参考文献[36]。 -
ERTBP考虑了主天体公转轨道的偏心率，相比CRTBP可以更精确地描述三体系统内航天器的运动情况。当两主天体围绕公共质心作椭圆轨道运动时，二者之间的距离呈现周期性变化，不再为一固定值，此时仍以时间作为独立变量不再适合，转而以主天体公转轨道的真近点角
$f$ 作为独立变量。相应的质量$[M]$ 、长度$[L]$ 和时间$[T]$ 的归一化单位取为$$ \left\{\begin{split} & {[M]=m_{1}+m_{2}} \\ & {[L]=L_{12}=\dfrac{a\left(1-e^{2}\right)}{1+e \cos f}} \\ & {[T]=\sqrt{\dfrac{L_{12}^{3}}{G\left(m_{1}+m_{2}\right)}}=\dfrac{\sqrt{1+e \cos f}}{\dot{f}}} \end{split}\right. $$ (4) 其中：
$\hat R$ 为当前时刻主天体之间的瞬时距离；$a$ 为主天体公转轨道的半长轴；$e$ 为公转轨道的偏心率；$f$ 为公转轨道的真近点角；$\dot f$ 为公转轨道的真近点角变化率。在ERTBP中，为研究和计算的方便，通常选择在脉动坐标系下进行研究。脉动坐标系的原点同样位于两主天体的质心，三轴指向与会合坐标系的三轴指向相同，但由于此坐标系的旋转角速率是变化的，且单位长度对应的真实物理长度也是时变的，为体现该坐标系的这种脉动特性，故称为脉动坐标系[19]。因此，在脉动坐标系下的航天器运动的动力学模型为
$$\left\{ \begin{split} & \ddot x - 2\dot y = \dfrac{1}{{1 + e\cos f}}{\varOmega _x} \\ &\ddot y + 2\dot x = \dfrac{1}{{1 + e\cos f}}{\varOmega _y} \\ &\ddot z + z = \dfrac{1}{{1 + e\cos f}}{\varOmega _z} \\ \end{split} \right.$$ (5) 其中：
$\varOmega $ 为ERTBP下的有效势函数，定义如下$$\varOmega = \dfrac{1}{2}({x^2} + {y^2} + {z^2}) + \dfrac{{1 - \mu }}{{{r_1}}} + \dfrac{\mu }{{{r_2}}}$$ (6) 而
${\varOmega _x},{\varOmega _y},{\varOmega _z}$ 为势函数$\varOmega $ 关于三轴的偏导数，具体表达式为$$\left\{ \begin{split} & {\varOmega _x} = x - (1 - \mu )\dfrac{{x + \mu }}{{r_1^3}} - \mu \dfrac{{x + \mu - 1}}{{r_2^3}} \\ & {\varOmega _y} = y - (1 - \mu )\dfrac{y}{{r_1^3}} - \mu \dfrac{y}{{r_2^3}} \\ & {\varOmega _z} = z - (1 - \mu )\dfrac{z}{{r_1^3}} - \mu \dfrac{z}{{r_2^3}} \\ \end{split} \right.$$ (7) 将上述偏导数代入动力学方程（5）可得
$$ \left\{ {\begin{split} & {\ddot x - 2\dot y = \dfrac{1}{{1 + e\cos f}}\left[ {x - (1 - \mu )\dfrac{{x + \mu }}{{r_1^3}} - \mu \dfrac{{x + \mu - 1}}{{r_2^3}}} \right]}\\ & {\ddot y + 2\dot x = \dfrac{1}{{1 + e\cos f}}\left[ {y - (1 - \mu )\dfrac{y}{{r_1^3}} - \mu \dfrac{y}{{r_2^3}}} \right]}\\ & {\ddot z + z = \dfrac{1}{{1 + e\cos f}}\left[ {z - (1 - \mu )\dfrac{z}{{r_1^3}} - \mu \dfrac{z}{{r_2^3}}} \right]} \end{split}} \right. $$ (8) 为获得ERTBP下平动点附近的航天器运动近似解析解，首先需要计算ERTBP下平动点的位置。由于平动点为三体系统的动平衡点，在平动点处满足速度、加速度为零，即
$$\dot x = \dot y = \dot z = \ddot x = \ddot y = \ddot z = 0$$ (9) 将此条件代入式（8），可得
$$ \left\{ {\begin{split} & {\dfrac{1}{{1 + e\cos f}}\left[ {x - (1 - \mu )\dfrac{{x + \mu }}{{r_1^3}} - \mu \dfrac{{x + \mu - 1}}{{r_2^3}}} \right] = 0}\\ & {\dfrac{1}{{1 + e\cos f}}\left[ {y - (1 - \mu )\dfrac{y}{{r_1^3}} - \mu \dfrac{y}{{r_2^3}}} \right] = 0}\\ & {\dfrac{1}{{1 + e\cos f}}\left[ {z - (1 - \mu )\dfrac{z}{{r_1^3}} - \mu \dfrac{z}{{r_2^3}}} \right] - z = 0} \end{split}} \right. $$ (10) 由于在Z轴上存在如下约束关系
$$\left\{ \begin{split} & z(e\cos f + \dfrac{{1 - \mu }}{{r_1^3}} + \dfrac{\mu }{{r_2^3}}) = 0 \\ & e\cos f + \dfrac{{1 - \mu }}{{r_1^3}} + \dfrac{\mu }{{r_2^3}} \ne 0 \\ \end{split} \right.$$ (11) 因此必有
$$z = 0$$ (12) 故ERTBP中的平动点位置和CRTBP一样，均在XY平面内，更具体的有
$$\left\{ \begin{split} & x - (1 - \mu )\dfrac{{x + \mu }}{{r_1^3}} - \mu \dfrac{{x + \mu - 1}}{{r_2^3}} = 0 \\ & y - (1 - \mu )\dfrac{y}{{r_1^3}} - \mu \dfrac{y}{{r_2^3}} = 0 \\ \end{split} \right.$$ (13) 这与CRTBP下平动点的计算公式完全相同。据此可知，在每一瞬时，ERTBP中的真近点角
$f$ 为一定值时，在脉动坐标系下的ERTBP平动点位置与会合坐标系下的CRTBP平动点位置完全相同，这也与文献[16]中的论述相符合，具体的平动点位置计算公式可以参考文献[36]，本文不再赘述。为获得ERTBP下平动点附近运动的近似解析解，需要在平动点
${x_{Li}}$ 处将ERTBP下的航天器运动动力学模型进行线性化$$ \delta \dot {{x}} = {\left. {\dfrac{{\partial f({{x}})}}{{\partial {{x}}}}} \right|_{{{x}} = {{{x}}_{{L_i}}}}}\delta {{x}} = {D_{{L_i}}}\delta {{x}} + {\rm H.O.T} $$ (14) 其中，H.O.T(High-Order-Terms)表示线性化过程中的高阶项。将线性化动力学模型式（14）按照三轴位置、速度展开可得式（15）。其中，
${\varOmega _{i,j}}(i,j = x,y,z)$ 是势函数$\varOmega $ 关于三轴的二阶偏导数，具体表达式为$$\delta \dot x = \left[ {\begin{array}{cccccccc} 0 & 0 &0&1&0&0 \\ 0 & 0 &0&0&1&0 \\ 0& 0&0&0&0&1 \\ {{\varOmega _{x,x}}}&{{\varOmega _{x,y}}}&{{\varOmega _{x,z}}}&0&2&0 \\ {{\varOmega _{y,x}}}&{{\varOmega _{y,y}}}&{{\varOmega _{y,z}}}&{ - 2}&0&0 \\ {{\varOmega _{z,x}}}&{{\varOmega _{z,y}}}&{{\varOmega _{z,z}}}&0&0&0 \end{array}} \right]\delta x$$ (15) $$\left\{ \begin{array}{l} {\varOmega _{x,x}} = \dfrac{1}{{1 + e\cos f}}\Bigg[1 - \dfrac{{1 - \mu }}{{r_1^3}} + 3\dfrac{{(1 - \mu ){{(x + \mu )}^2}}}{{r_1^5}} - \dfrac{\mu }{{r_2^3}} + 3\dfrac{{\mu {{(x + \mu - 1)}^2}}}{{r_2^5}}\Bigg] \\ {\varOmega _{y,y}} = \dfrac{1}{{1 + e\cos f}}\Bigg[1 - \dfrac{{1 - \mu }}{{r_1^3}} + 3\dfrac{{(1 - \mu ){y^2}}}{{r_1^5}} - \dfrac{\mu }{{r_2^3}} + 3\dfrac{{\mu {y^2}}}{{r_2^5}}\Bigg] \\ {\varOmega _{z,z}} = \dfrac{1}{{1 + e\cos f}}\Bigg[ - \dfrac{{1 - \mu }}{{r_1^3}} + 3\dfrac{{(1 - \mu ){z^2}}}{{r_1^5}} - \dfrac{\mu }{{r_2^3}} + 3\dfrac{{\mu {z^2}}}{{r_2^5}}\Bigg] - \dfrac{{e\cos f}}{{1 + e\cos f}} \\ {\varOmega _{x,y}} = {\varOmega _{y,x}} = \dfrac{1}{{1 + e\cos f}}\Bigg[3\dfrac{{(1 - \mu )(x + \mu )y}}{{r_1^5}} + 3\dfrac{{\mu (x + \mu - 1)y}}{{r_2^5}}\Bigg] \\ {\varOmega _{x,z}} = {\varOmega _{z,x}} = \dfrac{1}{{1 + e\cos f}}\Bigg[3\dfrac{{(1 - \mu )(x + \mu )z}}{{r_1^5}} + 3\dfrac{{\mu (x + \mu - 1)z}}{{r_2^5}}\Bigg] \\ {\varOmega _{y,z}} = {\varOmega _{z,y}} = \dfrac{1}{{1 + e\cos f}}\Bigg[3\dfrac{{(1 - \mu )yz}}{{r_1^5}} + 3\dfrac{{\mu yz}}{{r_2^5}}\Bigg] \\ \end{array} \right.$$ (16) 对于共线平动点，满足
$y = z = 0$ ，将此条件代入上述线性化动力学模型式（15），化简后可得$$\delta \dot x = \left[ {\begin{array}{cccccccc} 0&0&0&1&0&0 \\ 0&0&0&0&1&0 \\ 0&0&0&0&0&1 \\ {g(1 + 2\bar \mu )}&0&0&0&2&0 \\ 0&{g(1 - \bar \mu )}&0&{ - 2}&0&0 \\ 0&0&{g( - \bar \mu ) + k}&0&0&0 \end{array}} \right]\delta x$$ (17) 其中：参数
$g$ ，$\bar \mu $ ，$k$ 的表达式如式（18）所示$$\left\{ \begin{array}{l} g = \dfrac{1}{{1 + e\cos f}} \\ \bar \mu = \dfrac{{1 - \mu }}{{{{\left| {x + \mu } \right|}^3}}} + \dfrac{\mu }{{{{\left| {x + \mu - 1} \right|}^3}}} > 1 \\ k = \dfrac{{ - e\cos f}}{{1 + e\cos f}} \\ \end{array} \right.$$ (18) 通过线性系统理论方法，计算得到线性化常微分方程组式（15）的6个特征值
${\lambda _i}(i = 1,\cdots,6)$ 为$$\left\{ \begin{array}{l} {\lambda _1} = \sqrt {\dfrac{{2g + g\bar \mu - 4 + \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} }}{2}} \\ {\lambda _2} = - {\lambda _1} \\ {\lambda _3} = i\sqrt {\dfrac{{ - 2g - g\bar \mu + 4 + \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} }}{2}} \\ {\lambda _4} = - {\lambda _3} \\ {\lambda _5} = i\sqrt {g\bar u - k} \\ {\lambda _6} = - {\lambda _5} \\ \end{array} \right.$$ (19) 对应的6个特征向量
${\vec v_i}(i = 1,\cdots,6)$ 为$$\left\{ \begin{array}{l} {{ v}_1} = \left[ {\begin{array}{*{20}{c}} 1&{{\eta _1}}&0&{{\lambda _1}}&{{\eta _2}}&0 \end{array}} \right] \\ {{ v}_2} = \left[ {\begin{array}{*{20}{c}} 1&{ - {\eta _1}}&0&{ - {\lambda _1}}&{{\eta _2}}&0 \end{array}} \right] \\ {{ v}_3} = \left[ {\begin{array}{*{20}{c}} 1&{{\eta _3}}&0&{{\lambda _3}}&{{\eta _4}}&0 \end{array}} \right] \\ {{ v}_4} = \left[ {\begin{array}{*{20}{c}} 1&{ - {\eta _3}}&0&{ - {\lambda _3}}&{{\eta _4}}&0 \end{array}} \right] \\ {{ v}_5} = \left[ {\begin{array}{*{20}{c}} 0&0&1&0&0&{{\lambda _5}} \end{array}} \right] \\ {{ v}_6} = \left[ {\begin{array}{*{20}{c}} 0&0&1&0&0&{ - {\lambda _5}} \end{array}} \right] \\ \end{array} \right.$$ (20) 其中，
${\eta _1},{\eta _2},{\eta _3},{\eta _4}$ 的计算公式如式（21）所示$$\left\{ \begin{array}{l} {\eta _1} = \dfrac{{\sqrt 2 \sqrt {2g + g\bar u + \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} - 4} \cdot [2g + g\bar u - \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} - 4]}}{{g(\bar u - 1)(3g\bar u + \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} + 4)}}\\ {\eta _2} = - \dfrac{{{\rm{4}}g{\rm{(2}}\bar u{\rm{ + 1)}}}}{{3g\bar u + \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} + 4}}\\ {\eta _3} = \dfrac{{\sqrt 2 \sqrt {2g + g\bar u - \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} - 4} \cdot (2g + g\bar u + \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} - 4)}}{{g(\bar u - 1)(3g\bar u - \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} + 4)}}\\ {\eta _4} = - \dfrac{{{\rm{4}}g{\rm{(2}}\bar u{\rm{ + 1)}}}}{{3g\bar u + \sqrt {9{g^2}{{\bar u}^2} - 8g\bar u - 16g + 16} + 4}} \end{array} \right.$$ (21) 由于上述6个特征向量线性无关，由此得到了线性化常微分方程组式（15）的6个线性无关的解，结合6个特征值，构成了线性化动力学模型的基本解集，方程组的通解可表示为
$$\left\{ \begin{array}{l} \delta x(t) = {A_1}{e^{{\lambda _1}t}} + {A_2}{e^{ - {\lambda _1}t}} + {A_3}\cos (\rm{Im} ({\lambda _3})t) +\\ \quad\quad\quad {A_4}\sin (\rm{Im} ({\lambda _3})t) \\ \delta y(t) = {k_1}{A_1}{e^{{\lambda _1}t}} - {k_1}{A_2}{e^{ - {\lambda _1}t}} - {k_2}{A_3}\sin (\rm{Im} ({\lambda _3})t)+ \\ \quad\quad\quad {k_2}{A_4}\cos (\rm{Im} ({\lambda _3})t) \\ \delta z(t) = {A_5}\cos (\rm{Im} ({\lambda _5})t) + {A_6}\sin (\rm{Im} ({\lambda _5})t) \\ \end{array} \right.$$ (22) 其中，
${A_i}(i = 1, \cdots,6)$ 同样是由初始条件确定的积分常数，参数${k_1}$ ，${k_2}$ 的取值为$${k_1} = {m_1},{k_2} = \rm{Im} ({m_3})$$ (23) 式（22）即为本文研究的ERTBP下共线平动点附近运动的近似解析解。
通过对比发现ERTBP下平动点附近运动的近似解析解与CRTBP下的近似解析解形式完全相同，二者的区别仅仅是特征值
${\lambda _1},{\lambda _3},{\lambda _5}$ 和相关系数${k_1},{k_2}$ 的计算公式不同而已，这也在某种程度上体现了CRTBP与ERTBP的一致性。 -
为验证所推导的ERTBP下的共线平动点附近运动近似解析解的正确性和精度，本文将推导的近似解析解式（22）与经典的CRTBP下的近似解析解式（3）进行了仿真对比。
首先对解析解的正确性进行验证。具体过程如下所述。由于ERTBP是更为一般的三体模型，可以在一定条件下退化为CRTBP这一特殊情况，因此只需要将式（22）中的主天体公转轨道的偏心率
$e$ 设置为0，即$e = 0$ ，这样便将ERTBP退化为CRTBP，然后将退化的ERTBP平动点附近运动近似解析解式（22）与经典的CRTBP平动点附近运动近似解析解式（3）进行同样参数[A1，A2，A3，A4，A5，A6]下的仿真对比，便可以验证本文推导的近似解析解的正确性。本文以地–月系统L2平动点附近的Lissajous轨道作为仿真对象，在归一化单位DU下，相应的参数[A1，A2，A3，A4，A5，A6]设置为[0，0，0.05，0.05，0.01，0.01]，仿真时间设置为 [0，10]。仿真结果如图2和图3所示。
图 2ERTBP下的近似解析解与CRTBP下的近似解析解生成的Lissajous轨道对比
Figure 2.The comparison of Lissajous orbits generated by solutions in ERTBP and CRTBP
图 3ERTBP下的近似解析解相对CRTBP下的近似解析解的精度
Figure 3.The precision of the solutions in ERTBP relative to the solutions in CRTBP
根据仿真结果可知，退化后的ERTBP共线平动点附近运动近似解析解与经典的CRTBP共线平动点附近运动近似解析解高度一致，体现在Lissajous轨道的仿真结果上，本文推导的近似解析解的三轴最大绝对位置误差为[0，1.110 2 × 10-16，0]，最大相对位置误差为[0，3.979 6 × 10-14，0]，最大绝对速度误差为[0，1.110 2 × 10-16，0]，最大相对速度误差为[0，8.183 7 × 10-14，0]。仿真结果证明了本文推导的ERTBP共线平动点附近运动近似解析解的正确性。
为进一步分析ERTBP解析解与CRTBP解析解的精度差别，需要将主天体公转轨道的偏心率
$e$ 考虑进来。对于地月系统，真实偏心率$e = 0.054 \;88$ ，但为了定性地突出两解析解之间的精度差别，将偏心率$e$ 人为设置为更大的值$e = 0.3$ ，参数[A1，A2，A3，A4，A5，A6]同样设置为[0，0，0.05，0.05，0.01，0.01]，仿真时间设置为 [0,0.5]。仿真结果如图4和图5所示，由于篇幅限制，这里仅给出XY平面内的运动情况。
图 4CRTBP解析解精度
Figure 4.The precision of CRTBP analytical solutions
图 5ERTBP解析解精度
Figure 5.The precision of ERTBP analytical solutions
从图4和图5可以看出，分别利用CRTBP解析解和ERTBP解析解提供的初值，对非线性动力学模型（5）进行数值积分，均会出现发散现象，这是由三体问题的本质所决定的。然而，ERTBP解析解提供的初值经过数值积分后发散趋势更小，数值积分生成的“真实”轨道可以更好地维持在解析解生成的“理论”轨道附近。仿真结果证明了ERTBP解析解提供的初值相比CRTBP解析解提供的初值具有更高的精度。
-
为了更好地为未来深空探测任务服务和进一步完善椭圆型限制性三体问题下平动点附近运动的解析研究，本文参考借鉴圆型限制性三体问题下共线平动点附近运动近似解析解的研究方法，首先根据平动点的动平衡性质获得了平动点的位置信息，然后将椭圆型限制性三体问题下的非线性动力学模型在共线平动点处线性化，最后结合线性系统理论求解线性常微分方程组，获得了比圆型限制性三体问题更具一般性的椭圆型限制性三体问题下共线平动点附近运动的近似解析解。最后以地–月三体系统L2平动点Lissajous轨道为例，与经典的圆型限制性三体问题下的近似解析解进行仿真对比，结果证明了本文推导的椭圆型限制性三体问题下的共线平动点附近运动的近似解析解的正确性，同时也表明椭圆型限制性三体问题近似解析解相比圆型限制性三体问题近似解析解具有更高的精度。
Approximate Analytical Solutions of Motion near the Collinear Libration-Points in Restricted Three-Body Problem
-
摘要:随着深空探测成为航天领域的研究热点，与其密切相关的三体问题基础研究也日益重要，尤其是在深空探测任务设计中处于基础地位的共线平动点附近运动的研究，更是具有重要的工程应用价值。在圆型限制性三体问题下，对共线平动点附近运动近似解析解的研究已经较为全面，但在更接近真实情况、更具一般性的椭圆型限制性三体问题下，相应的研究却相对较少。针对此背景，参考借鉴圆型限制性三体问题的研究方法，首先根据平动点的特性计算出平动点的位置，然后将非线性三体动力学模型在共线平动点处线性化，最后结合线性系统理论，获得了椭圆型限制性三体问题下共线平动点附近运动的近似解析解，并将其与经典的圆型限制性三体问题下的近似解析解进行对比分析，仿真结果证明了方法的有效性，同时也表明所推导的椭圆型限制性三体问题解析解相比圆型限制性三体问题解析解具有更高的精度。
-
关键词:
- 圆型限制性三体问题/
- 椭圆型限制性三体问题/
- 非线性系统/
- 共线平动点/
- 解析解
Abstract:With deep space exploration becoming a research focus in aerospace, corresponding fundamental research on three-body problem is of increasingly significant, especially the motion analysis near the collinear libration-points, which play a leading role in deep space mission design The approximate analytical solutions of motion near the collinear libration-points in circular restricted three-body problem has been obtained, however, there are relatively fewer studies about the solutions in elliptic restricted three-body problem, although it is more realistic and general than circular restricted three-body problem. Based on it, the approximate analytical solutions of motion near the collinear libration-points in elliptic restricted three-body problem are deduced by referencing the method used in circular restricted three-body problem, the positions of libration-points are obtained according to its characteristic, then the nonlinear dynamic model is linearized at the collinear libration-points, and the approximate analytical solutions are finally obtained using the linear system theory and compared with the solutions of the circular restricted three-body problem. Simulation results indicated the method is valid and the deduced analytical solutions have higher precision than that of the circular restricted three-body problem.Highlights● The positions of libration points in elliptic restricted three-body problem are given. ● The approximate analytical solutions of motion around collinear pibration points in elliptic restricted three-body problem are derived. ● The solutions for elliptic restricted three-body problem have higher precision than that of the circular restricted three-body problem was proved. -
图 1质心旋转坐标系
$O - XYZ$ 和平动点旋转坐标系${L_2} - xyz$ Fig. 1Barycentric rotating coordinate system
$O - XYZ$ and libration-point centered rotating coordinate system${L_2} - xyz$
图 2ERTBP下的近似解析解与CRTBP下的近似解析解生成的Lissajous轨道对比
Fig. 2The comparison of Lissajous orbits generated by solutions in ERTBP and CRTBP
图 3ERTBP下的近似解析解相对CRTBP下的近似解析解的精度
Fig. 3The precision of the solutions in ERTBP relative to the solutions in CRTBP
-
[1] 孟云鹤, 张跃东, 陈琪锋. 平动点航天器动力学与控制[M]. 北京: 科学出版社, 2015. [2] 乔栋,崔平远,尚海滨. 基于椭圆型限制性三体模型的借力飞行机理研究[J]. 宇航学报,2010,31(1):36-43.doi:10.3873/j.issn.1000-1328.2010.01.005QIAO D,CUI P Y,SHANG H B. Research on gravity-assist mechanism in elliptic restricted Three-body Model[J]. Journal of Astronautics,2010,31(1):36-43.doi:10.3873/j.issn.1000-1328.2010.01.005 [3] GOMEZ G. Dynamics and mission design near libration points. Vol. 1. fundamentals: the case of collinear libration points （world scientific monograph series in mathematics; Vol. 2）[M]. Singapore: World Scientific, 2001. [4] RICHARDSON D L. Analytic construction of periodic orbits about the collinear points[J]. Celestial Mechanics,1980,22:241-253.doi:10.1007/BF01229511 [5] MASDEMONT J. High-order expansions of invariant manifolds of libration point orbits with application to mission design[J]. Dynamical Systems,2005,20(1):59-113.doi:10.1080/14689360412331304291 [6] JORBA À,MASDEMONT J. Dynamics in the center manifold of the collinear points of the restricted three body problem[J]. Physica. Section D:Nonlinear Phenomena,1999,132(1-2):189-213.doi:10.1016/S0167-2789(99)00042-1 [7] RITHWIK N, RAMANAN R V. Halo orbit design around Lagrangian points using gradient and non-gradient based optimization techniques[C]//AIP Conference Proceedings. [S.l.]: AIP Publishing, 2018, 1975（1）: 030019. [8] 刘刚,黄静,郭思岩. 星历模型地月系统平动点拟周期轨道设计研究[J]. 上海航天,2017(5):16-22.LIU G,HUANG J,GUO S Y. Quasi-periodic orbit design around Earth-Moon libration points using ephemeris model[J]. Aerospace Shanghai,2017(5):16-22. [9] 陶雪峰,连一君. 基于参数优化的平动点轨道族简单生成方法[J]. 中国科学:技术科学,2017(1):93-106.TAO X F,LIAN Y J. A simple method for generating families of libration points orbits based on parameters optimization[J]. Science Sinica Technologica,2017(1):93-106. [10] LARA M,PEREZ I L,LOPEZ R. Higher order approximation to the Hill problem dynamics about the libration points[J]. Communications in Nonlinear Science and Numerical Simulation,2018,59:612-628.doi:10.1016/j.cnsns.2017.12.007 [11] 郑越,泮斌峰,唐硕. 一种计算三体问题周期轨道的新方法[J]. 宇航学报,2017,38(4):384-392.ZHENG Y,PAN B F,TANG S. A novel method of periodic orbit computation for three-body problem[J]. Journal of Astronautics,2017,38(4):384-392. [12] LEI H L,XU B. High-order analytical solutions around triangular libration points in circular restricted three-body problem[J]. Monthly Notices of the royal Astronomical Society,2013,434(2):1376-1386.doi:10.1093/mnras/stt1099 [13] LIANG Y,XU M,XU S. High-order solutions of motion near triangular libration points for arbitrary value of μ[J]. Nonlinear Dynamics,2018,93(1):1-24.doi:10.1007/s11071-018-4315-x [14] QIAN Y J,YANG X D,ZHAI G Q,et al. Analytical and numerical construction of vertical periodic orbits about triangular libration points based on polynomial expansion relations among directions[J]. Astrophysics and Space Science,2017,362(8):136.doi:10.1007/s10509-017-3115-y [15] 翟冠峤,张伟,钱霙婧. 基于多项式展开的三角平动点垂直周期轨道解析构建方法研究[J]. 动力学与控制学报,2018,16(1):26-34.doi:10.6052/1672-6553-2017-52ZHAI G Q,ZHANG W,QIAN Y J. Construction methods of vertical periodic orbit around the triangular libration points based on polynomial expansion[J]. Journal of Dynamics and Control,2018,16(1):26-34.doi:10.6052/1672-6553-2017-52 [16] 刘林, 侯锡云. 深空探测器轨道力学[M]. 北京: 电气工业出版社, 2015. [17] 邓辉. 地月系共线平动点动力学特征与应用研究[D]. 南京: 南京大学, 2017.DENG H. Research on the dynamics and applications of the collinear libration point in the Earth-Moon system[D]. Nanjing: Nanjing University, 2017. [18] 熊欢欢,高有涛. 椭圆限制性三体问题模型下平动点拟周期轨道卫星的自主定轨分析[J]. 中国科技论文,2015,10(4):478-487.doi:10.3969/j.issn.2095-2783.2015.04.023XIONG H H,GAO Y T. Study of satellite autonomous orbit determination in libration point quasi-periodic orbit based on elliptic restricted three-body problem model[J]. China Sciencepaper,2015,10(4):478-487.doi:10.3969/j.issn.2095-2783.2015.04.023 [19] 卢松涛,赵育善. 太阳辐射对ER3BP拉格朗日点的影响[J]. 飞行力学,2009,27(5):70-74.LU S T,ZHAO Y S. Effect of solar radiation pressure on the ER3BP Lagrange points[J]. Flight Dynamics,2009,27(5):70-74. [20] 雷汉伦. 平动点、不变流形及低能轨道[D]. 南京: 南京大学, 2015.LEI H L. Equilibrium point, Invariant manifold and low-energy trajectory[D]. Nanjing: Nanjing University, 2015. [21] HOU X Y,LIU L. On motions around the collinear libration points in the elliptic restricted three-body problem[J]. Monthly Notices of the Royal Astronomical Society,2011,415(4):3552-3560.doi:10.1111/j.1365-2966.2011.18970.x [22] LEI H,XU B,HOU X,et al. High-order solutions of invariant manifolds associated with libration point orbits in the elliptic restricted three-body system[J]. Celestial Mechanics and Dynamical Astronomy,2013,117(4):349-384.doi:10.1007/s10569-013-9515-6 [23] HOU X Y,LIU L. On quasi-periodic motions around the triangular libration points of the real Earth-Moon system[J]. Celestial Mechanics and Dynamical Astronomy,2010,108(3):301-313.doi:10.1007/s10569-010-9305-3 [24] LEI H,XU B. High-order solutions around triangular libration points in the elliptic restricted three-body problem and applications to low energy transfers[J]. Communications in Nonlinear Science and Numerical Simulation,2014,19(9):3374-3398.doi:10.1016/j.cnsns.2014.01.019 [25] GONG S P,LI J F. Solar sail periodic orbits in the elliptic restricted three-body problem[J]. Celestial Mechanics and Dynamical Astronomy,2015,121(2):121-137.doi:10.1007/s10569-014-9590-3 [26] KOH D,ANDERSON R L. Periodic orbit-attitude solutions in the planar elliptic restricted three-body problem[J]. Journal of Guidance Control and Dynamics,2018,41(3):644-657.doi:10.2514/1.G002885 [27] GOMEZ G,MASDEMONTJ J,SIMO C. Quasi-halo orbits associated with libration points[J]. Journal of the Astronautical Sciences,1998,46(2):135-176. [28] NARAYAN A,SINGH N. Characteristics exponents of the triangular solutions in the elliptical restricted three-body problem under radiating primaries[J]. Differential Equations and Dynamical Systems,2016,24(3):329-343.doi:10.1007/s12591-014-0230-x [29] FARQUHAR R W,KAMEL A A. Quasi-periodic orbits about the translunar libration point[J]. Celestial Mechanics,1973,7(4):458-473.doi:10.1007/BF01227511 [30] SINGH J,GYEGWE J. Analytic approximation solutions of lyapunov orbits around the collinear equilibrium points for binary α-centuari system:the planar case[J]. British Journal of Mathematics and Computer Science,2017,22(1):1-18. [31] HOU X Y,LIU L. On quasi-periodic motions around the collinear libration points in the real Earth-Moon system[J]. Celestial Mechanics and Dynamical Astronomy,2011,110(1):71-98.doi:10.1007/s10569-011-9340-8 [32] ANSARI A A. Periodic orbits around lagrangian points of the circular restricted four-body problem[J]. Invertis Journal of Science and Technology,2014,7(1):29-38. [33] PAPADAKIS K E. Families of three-dimensional periodic solutions in the circular restricted four-body problem[J]. Astrophysics and Space Science,2016,361(4):1-14. [34] JAGADISH S,PERDIOU A E,MRUMUN G J,et al. Periodic solutions around the collinear equilibrium points in the perturbed restricted three-body problem with triaxial and radiating primaries for binary HD 191408,Kruger 60 and HD 155876 systems[J]. Applied Mathematics and Computation,2018,325:358-374.doi:10.1016/j.amc.2017.11.052 [35] HOWELL K C,PERNICKA H J. Numerical determination of Lissajous trajectories in the restricted three-body problem[J]. Celestial mechanics,1987,41(1-4):107-124.doi:10.1007/BF01238756 [36] 赵育善, 师鹏, 张晨.深空飞行动力学[M]. 北京: 中国宇航出版社, 2016. -
点击查看大图
图(5)
计量
- 文章访问数:1031
- HTML全文浏览量:264
- PDF下载量:165
- 被引次数:0
