CN103488830A - 一种基于Cycler轨道的地月往返的任务仿真*** - Google Patents

Abstract

本发明公开了一种基于Cycler轨道的地月往返的任务仿真***，该***能够对地球与月球之间往返的任务进行仿真，横跨了航天器轨道设计，航天器任务分析以及航天器动力学仿真等。本发明通过构建共振Cycler轨道模型、构建同宿Cycler轨道模型、多重打靶法轨道修正模块和Lambert转移轨道获取模块来对共振型Cycler轨道及同宿型Cycle轨道设计修正、Lambert转移轨道设计修正，发射交会窗口分析以及地月往返***任务分析的技术问题。本发明***在获得地月Cycler轨道时考虑了太阳引力摄动，并且运用多重打靶法得到近似周期Cycler轨道，与现有方法计算的周期Cycler轨道相比精度有所提高。

Description

一种基于Cycler轨道的地月往返的任务仿真***

技术领域

本发明涉及一种地球与月球之间往返的任务仿真，更特别地说，是指一种基于Cycler轨道的地月往返的任务仿真***。

背景技术

中国首次月球探测工程由月球探测卫星、运载火箭、发射场、测控和地面应用等五大***组成，中国探月卫星工程还有五大工程目标：一是研制和发射中国第一颗探月卫星；二是初步掌握绕月探测基本技术；三是首次开展月球科学探测；四是初步构建月球探测航天工程***；五是为月球探测后续工程积累经验。为此要突破月球探测卫星的关键技术；初步建立中国的深空探测工程大***；验证有效载荷和数据解译等各项关键技术；初步建立中国深空探测技术研制体系；培养相应的人才队伍。

中国首次月球探测工程四大科学任务为：

一、获取月球表面三维立体影像，精细划分月球表面的基本构造和地貌单元，进行月球表面撞击坑形态、大小、分布、密度等的研究，为类地行星表面年龄的划分和早期演化历史研究提供基本数据，并为月面软着陆区选址和月球基地位置优选提供基础资料等。

二、分析月球表面有用元素含量和物质类型的分布特点，主要是勘察月球表面有开发利用价值的钛、铁等14种元素的含量和分布，绘制各元素的全月球分布图，月球岩石、矿物和地质学专题图等，发现各元素在月表的富集区，评估月球矿产资源的开发利用前景等。

三、探测月壤厚度，即利用微波辐射技术，获取月球表面月壤的厚度数据，从而得到月球表面年龄及其分布，并在此基础上，估算核聚变发电燃料氦3的含量、资源分布及资源量等。

四、探测地球至月球的空间环境。月球与地球平均距离为38万公里，处于地球磁场空间的远磁尾区域，卫星在此区域可探测太阳宇宙线高能粒子和太阳风等离子体，研究太阳风和月球以及地球磁场磁尾与月球的相互作用。

Cycler轨道是指周期性往返于地球和月球之间，在行星附近绕飞而不停留的轨道。运行于循环轨道上的飞行器可长时间(几年甚至十几年)保持行星间飞行而不需要轨道机动(或者只需要很小的轨道机动)，因而被认为是节约能量的一种经济型长期任务方式。Cycler轨道方案按照绕日圈数不同还可再细分为整圈轨道、半圈轨道、一般轨道等。Cycler轨道可以分为共振型Cycler轨道及同宿型Cycler轨道。

圆型限制性三体问题(Circular Restricted Three-Body Problem，CR3BP)描述质量相对无限小的第三体在两个围绕其公共质心做圆周运动的主天体的引力作用下的运动。参考2010年5月，李明涛的《中国科学院研究生院博士学位论文》，第19页至第21页的相关内容。

双圆模型(Bi-Circalar Model，BCM)是研究一个无限小质量体在月亮(绕***质心运行在圆轨道上)、太阳和地球(围绕共同质心运行在圆轨道上)的***的万有引力作用下的运动规律的基本模型。

发明内容

本发明的目的是提出一种基于Cycler轨道的地月往返的任务仿真***，该***在圆形限制性三体模型下产生共振型Cycler轨道和同宿型Cycler轨道。由于月球轨道的偏心率及其他天体(如太阳、木星等)引力影响，使得圆形限制性三体模型下建立的轨道与真实情况不同，甚至由于摄动的影响，这些轨道无法保持恒定。因此，以在圆形限制性三体模型下建立的轨道为初值，在双圆模型下进行优化，利用多重打靶法来修正共振型周期轨道和同宿型周期轨道。修正后的共振型周期轨道和同宿型周期轨道能够定期与地球、月球相遇，产生了以地球或月球为目标的发射窗口。为完成完整的地月转移轨道，本发明仿真***还根据二体Lambert问题的经典解法，生成了地球起飞、月球着陆的“地球-周期转移轨道”；月球起飞、地球着陆的“周期转移轨道-月球”；并在四体模型下进行了修正。在本发明仿真***中计算并比较这两种转移轨道的燃料消耗量及飞行时间。

本发明是一种基于Cycler轨道的地月往返的任务仿真***，该仿真***包括有构建共振Cycler轨道模型10、构建同宿Cycler轨道模型20、多重打靶法轨道修正模块30和Lambert转移轨道获取模块40。

所述构建共振Cycler轨道模型10第一方面依据CR3BP模型构建M_CR3BP动力学模型；第二方面采用BCM模型对所述的M_CR3BP动力学模型进行优化，得到M_BCM动力学模型。

所述构建同宿Cycler轨道模型20依据BCM模型构建M_L动力学模型。

所述多重打靶法轨道修正模块30第一方面采用多重打靶法分别对M_BCM进行修正，得到修正后的共振Cycler轨道动力学模型DM_BCM；第二方面对DM_BCM采用四阶龙格库塔法获得共振Cycler轨道SM_BCM；第三方面采用多重打靶法分别对M_L进行修正，得到修正后的同宿Cycler轨道动力学模型DM_L；第四方面对DM_L采用四阶龙格库塔法获得同宿Cycler轨道SM_L。

所述Lambert转移轨道获取模块40第一方面采用高斯-全局变量复合算法对SM_BCM进行处理，得到第一套Lambert转移轨道 $S_{\mathrm{BCM}}=\{s_1^{\mathrm{BCM}},s_2^{\mathrm{BCM}},s_3^{\mathrm{BCM}},s_4^{\mathrm{BCM}}\};$ 第二方面采用微分修正算法对进行修正，得到第一套修正后的Lambert转移轨道 ${\mathrm{DS}}_{\mathrm{BCM}}={\{\mathrm{ds}}_1^{\mathrm{BCM}},{\mathrm{ds}}_2^{\mathrm{BCM}},{\mathrm{ds}}_3^{\mathrm{BCM}},{\mathrm{ds}}_4^{\mathrm{BCM}}\};$ 第三方面依据发射点纬度和白赤交角变化规律对SM_BCM进行发射窗口及交会窗口的分析，获取航天器发射和入轨的时机；第四方面采用高斯-全局变量复合算法对SM_L进行处理，得到第二套Lambert转移轨道

第五方面采用微分修正算法对

进行修正，得到第二套修正后的Lambert转移轨道

第六方面通过比较地月往返任务的燃料消耗量、总飞行时间能够获知SM_BCM和SM_L轨道各自的优势，为小推力航天器的深空探测的地月往返任务提供优化设计指标。

本发明仿真***的优点在于：

①本发明在圆形限制性三体模型下产生共振型Cycler轨道和同宿型Cycler轨道。是为了修正现有共振型Cycler轨道及同宿型Cycle轨道中没有考虑太阳引力摄动对轨道设计的影响。

②本发明***在获得地月Cycler轨道时考虑了太阳引力摄动，并且运用多重打靶法得到近似周期Cycler轨道，与现有方法计算的周期Cycler轨道相比精度有所提高。

③本发明***在分析航天器由地面起飞至Cycler轨道的交会窗口和航天器由月面起飞至共振Cycler轨道的交会窗口时，提供一种迭代算法，可以在考虑太阳引力摄动情况下得到转移窗口，与现有Lambert方法计算结果相比精度有所提高。

附图说明

图1是本发明基于Cycler轨道的地月往返的任务仿真***的结构框图。

图2是地月质心惯性坐标系O-XYZ和地月质心旋转坐标系O-xyz的坐标系示意图。

图2A是日地旋转坐标系O_S-X_SY_SZ_S的坐标系示意图。

图3是经本发明的多重打靶法修正后的共振Cycler轨道动力学模型DM_BCM图。

图3A是经本发明的多重打靶法修正后的同宿Cycler轨道动力学模型DM_L图。

图4是Lambert问题下航天动力学中的两点边界值获取的示意图。

图5是经本发明高斯-全局变量复合算法处理SM_BCM后的仿真结果图。

图5A是经微分修正算法修正后的第一套修正Lambert转移轨道的仿真结果图。

图5B是经微分修正算法修正后的第二套修正Lambert转移轨道的仿真结果图。

图6是西昌直接入轨的窗口示意图。

图6A是最省燃料霍曼转移所需速度增量的轨道示意图。

图6B是最费燃料霍曼转移所需速度增量的轨道示意图。

图6C是在地月质心惯性坐标系O-XYZ下发射轨迹示意图。

图6D是地月质心旋转坐标系O-xyz下发射轨迹示意图。

图7是月面起飞的西昌直接入轨的窗口示意图。

图7A是停泊轨道高度为100km的约束条件下的搜索结果示意图。

图7B是地理经度为[-170°，-135°]∪[-43°，-5°]∪[45°，108°]∪[150°，165°]的入轨窗口示意图。

图7C是二体Lambert问题下的同宿Lambert轨道仿真结果图。

图7D是经微分修正算法后的第二套修正Lambert转移轨道 ${\mathrm{DS}}_L={\{\mathrm{ds}}_1^L,{\mathrm{ds}}_2^L\}$ 的仿真结果。

具体实施方式

下面将结合附图对本发明做进一步的详细说明。

参见图1所示，本发明是一种基于Cycler轨道的地月往返的任务仿真***，该仿真***包括有构建共振Cycler轨道模型10、构建同宿Cycler轨道模型20、多重打靶法轨道修正模块30和Lambert转移轨道获取模块40。

所述构建共振Cycler轨道模型10第一方面依据CR3BP模型构建M_CR3BP动力学模型；第二方面采用BCM模型对所述的M_CR3BP动力学模型进行优化，得到M_BCM动力学模型。

所述构建同宿Cycler轨道模型20依据BCM模型构建M_L动力学模型。

所述多重打靶法轨道修正模块30第一方面采用多重打靶法分别对M_BCM进行修正，得到修正后的共振Cycler轨道动力学模型DM_BCM；第二方面对DM_BCM采用四阶龙格库塔法获得共振Cycler轨道SM_BCM；第三方面采用多重打靶法分别对M_L进行修正，得到修正后的同宿Cycler轨道动力学模型DM_L；第四方面对DM_L采用四阶龙格库塔法获得同宿Cycler轨道SM_L。

所述Lambert转移轨道获取模块40第一方面采用高斯-全局变量复合算法对SM_BCM进行处理，得到第一套Lambert转移轨道 $S_{\mathrm{BCM}}=\{s_1^{\mathrm{BCM}},s_2^{\mathrm{BCM}},s_3^{\mathrm{BCM}},s_4^{\mathrm{BCM}}\};$ 第二方面采用微分修正算法对

进行修正，得到第一套修正后的Lambert转移轨道 ${\mathrm{DS}}_{\mathrm{BCM}}={\{\mathrm{ds}}_1^{\mathrm{BCM}},{\mathrm{ds}}_2^{\mathrm{BCM}},{\mathrm{ds}}_3^{\mathrm{BCM}},{\mathrm{ds}}_4^{\mathrm{BCM}}\};$ 第三方面依据发射点纬度和白赤交角变化规律对SM_BCM进行发射窗口及交会窗口的分析，获取航天器发射和入轨的时机；第四方面采用高斯-全局变量复合算法对SM_L进行处理，得到第二套Lambert转移轨道

第五方面采用微分修正算法对

进行修正，得到第二套修正后的Lambert转移轨道

第六方面通过比较地月往返任务的燃料消耗量、总飞行时间能够获知SM_BCM和SM_L轨道各自的优势，为小推力航天器的深空探测的地月往返任务提供优化设计指标。

(一)构建共振Cycler轨道模型10

在本发明中，依据CR3BP模型得到共振Cycler轨道动力学模型M_CR3BP：

$M_{\mathrm{CR}3\mathrm{BP}}=\left\{\begin{array}{c}\frac{\∂U}{\∂x}={\mathrm{HA}}^x-2V{A}^y\\ \frac{\∂U}{\∂y}={\mathrm{HA}}^y+2V{A}^x\\ \frac{\∂U}{\∂z}={\mathrm{HA}}^z\end{array}\right.---\left(1\right)$

如图2所示，地月质心惯性坐标系O-XYZ和地月质心旋转坐标系O-xyz。其中O为地月质心，以地月***运动平面为XY坐标面。地月质心惯性坐标系中X为初始时刻地球指向月球的方向，Z方向为地月***的角速度方向，Y方向依据X和Z方向右手规则确定；地月质心旋转坐标系中的z方向与地月质心惯性坐标系的Z方向重合，x方向始终为地球指向月球的方向，y根据右手系法则建立。设地球为P₁，月球为P₂，航天器为P，航天器P在地月质心旋转坐标系O-xyz下的位置记为(x_P，y_P，z_P)；航天器P到地球P₁的距离记为R₁( $R_1={\[{(x_P-{\mathrm{\μ}}_1)}^2+{y_P}^2+{z_P}^2\]}^{\frac{1}{2}}$ )，航天器P到月球P₂的距离记为R₂( $R_2={\[{(x_P+1-{\mathrm{\μ}}_1)}^2+{y_P}^2+{z_P}^2\]}^{\frac{1}{2}}$ )，航天器P到地月质心O的距离记为R。

在地月质心旋转坐标系O-xyz下，公式(1)字母的物理意义为：

表示x方向上的偏导数；

U表示航天器P的势函数，且 $U=\frac{1}{2}({x_P}^2+{x_P}^2)+\frac{1-{\mathrm{\μ}}_1}{R_1}+\frac{{\mathrm{\μ}}_1}{R_2},$ μ₁表示月球与地球的质量比，一般取值为0.01215；

表示y方向上的偏导数；

表示z方向上的偏导数；

VA^x表示航天器P在x方向上的速度；

VA^y表示航天器P在y方向上的速度；

HA^x表示航天器P在x方向上的加速度；

HA^y表示航天器P在y方向上的加速度；

HA^z表示航天器P在z方向上的加速度。

在本发明中，根据BCM模型对M_CR3BP进行修正，得到BCM模型下动力学模型M_BCM：

$M_{\mathrm{BCM}}=\left\{\begin{array}{c}{\mathrm{HB}}^{X_S}={2\mathrm{VB}}^{Y_S}+{\mathrm{xs}}_P-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{xs}}_P+{\mathrm{\μ}}_2}{r_{\mathrm{PS}}^3}\\ {\mathrm{\μ}}_2\frac{{\mathrm{xs}}_P-1+{\mathrm{\μ}}_2}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{xs}}_P-x_{P_2}}{r_{\mathrm{PM}}^3}\\ {\mathrm{HB}}^{Y_S}=-2V{B}^{X_S}+{\mathrm{ys}}_P-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{ys}}_P}{r_{\mathrm{PS}}^3}-\\ {\mathrm{\μ}}_2\frac{{\mathrm{ys}}_P}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{ys}}_P-y_{P_2}}{r_{\mathrm{PM}}^3}\\ {\mathrm{HB}}^{Z_S}=-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{zs}}_P}{r_{\mathrm{PS}}^3}-{\mathrm{\μ}}_2\frac{{\mathrm{zs}}_P}{r_{\mathrm{PM}}^3}-m_M\frac{{\mathrm{zs}}_P}{r_{\mathrm{PM}}^3}\end{array}\right.---\left(2\right)$

在BCM模型下，使用了日地旋转坐标系P_S-X_SY_SZ_S，如图2A所示，太阳为P₃，航天器P在日地旋转坐标系O_S-X_SY_SZ_S下的位置记为(xs_P，ys_P，zs_P)；月球P₂在日地旋转坐标系O_S-X_SY_SZ_S下的位置记为r_PS表示航天器P与太阳的无量纲化距离，且 $r_{\mathrm{PS}}={[{({\mathrm{xs}}_P+{\mathrm{\μ}}_2)}^2+{{\mathrm{ys}}_P}^2+{{\mathrm{zs}}_P}^2]}^{\frac{1}{2}};$ r_PE表示航天器P与地球的无量纲化距离，且 $r_{\mathrm{PE}}={[{({\mathrm{xs}}_P-1+{\mathrm{\μ}}_2)}^2+{{\mathrm{ys}}_P}^2+{{\mathrm{zs}}_P}^2]}^{\frac{1}{2}};$ r_PM表示航天器P与月球的无量纲化距离，且 $r_{\mathrm{PM}}={[{({\mathrm{xs}}_P-x_{P_2})}^2+{({\mathrm{ys}}_P-y_{P_2})}^2+{{\mathrm{zs}}_P}^2]}^{\frac{1}{2}};$ μ₂表示地球与太阳的质量比，一般取值为0.000003003，m_M表示月球无量纲化质量。以日地质心O_S为原点，日地连线为X_S轴，指向地球为X_S轴的正方向。Y_S轴垂直于日地运动平面，根据右手定则，可以确定Z_S轴的方向。

在日地旋转坐标系O_S-X_SY_SZ_S下，公式(2)中字母的物理意义为：

表示所在M_BCM模型下的航天器P在X_S方向上的速度；

表示所在M_BCM模型下的航天器P在Y_S方向上的速度；

表示所在M_BCM模型下的航天器P在X_S方向上的加速度；

表示所在M_BCM模型下的航天器P在Y_S方向上的加速度；

表示所在M_BCM模型下的航天器P在Z_S方向上的加速度。

(二)构建同宿Cycler轨道模型20

在本发明中，依据BCM模型得到同宿Cycler轨道动力学模型M_L：

$M_{\mathrm{BCM}}=\left\{\begin{array}{c}{\mathrm{HB}}^{X_S}={2\mathrm{VB}}^{Y_S}+{\mathrm{xs}}_P-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{xs}}_P+{\mathrm{\μ}}_2}{r_{\mathrm{PS}}^3}\\ {\mathrm{\μ}}_2\frac{{\mathrm{xs}}_P-1+{\mathrm{\μ}}_2}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{xs}}_P-x_{P_2}}{r_{\mathrm{PM}}^3}\\ {\mathrm{HB}}^{Y_S}=-2V{B}^{X_S}+{\mathrm{ys}}_P-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{ys}}_P}{r_{\mathrm{PS}}^3}-\\ {\mathrm{\μ}}_2\frac{{\mathrm{ys}}_P}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{ys}}_P-y_{P_2}}{r_{\mathrm{PM}}^3}\\ {\mathrm{HB}}^{Z_S}=-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{zs}}_P}{r_{\mathrm{PS}}^3}-{\mathrm{\μ}}_2\frac{{\mathrm{zs}}_P}{r_{\mathrm{PM}}^3}-m_M\frac{{\mathrm{zs}}_P}{r_{\mathrm{PM}}^3}\end{array}\right.---\left(3\right)$

在日地旋转坐标系O_S-X_SY_SZ_S下，公式(3)中字母的物理意义为：

表示所在M_L模型下的航天器P在X_s方向上的速度；

表示所在M_L模型下的航天器P在Y_S方向上的速度；

表示所在M_L模型下的航天器P在X_s方向上的加速度；

表示所在M_L模型下的航天器P在Y_S方向上的加速度；

表示所在M_L模型下的航天器P在Z_S方向上的加速度。

(三)多重打靶法轨道修正模块30

在本发明中，采用多重打靶法对初始化的M_BCM模型进行修正，得到修正后的共振Cycler轨道动力学模型DM_BCM。所述DM_BCM是以图形(如图3所示)进行表征的。所述多重打靶法参考1988年发表的《一类最佳飞行轨迹的计算》，第九卷，第一期，A21至A22页，作者王培德等。

在一个周期的M_BCM模型上取500个等时间间隔的样点，任意一样点记为ini_b(i)(i＝1，…，500)；根据M_BCM模型采用四阶龙格库塔法积分，得到任意一点的位置速度状态量记为ini_a(i)(i＝1，…，500)；ini_a(i，：)是一个500×6的矩阵。

在每个时间间隔上，选取100个样点，并对样点采用ODE45积分，取出末点的位置速度状态量

记末点的位置速度状态量

与任意一点的位置速度状态量ini_a(i，：)的差值为F(i，：)，即

应用差值F(i，：)是使得解算值和实际计算值能够无误差，在仿真过程中应尽量将误差缩减到最小。

对于上述差值F(i，：)采用牛顿迭代方法求极小值。所述牛顿迭代法参考2011年发表的《牛顿迭代法的MATLAB实现》，第六期，第20页，作者云磊。在设定精度为1×10^-10条件下，经过8次迭代，达到设定精度。

第1次迭代	2.4925020681966167×10-1
		第2次迭代	9.5043845337684230×10-2
第3次迭代	2.3381313323453772×10-2
		第4次迭代	1.9075592366893615×10-2
第5次迭代	4.5915032506800547×10-4
		第6次迭代	8.9114511750870658×10-5
第7次迭代	4.5782899616269825×10-9
		第8次迭代	8.9656802854146783×10-10

若大于所设定的精度，继续以牛顿迭代法重复上述迭代，直至小于所设定的精度为止。由此，可以得到修正后的模型DM_BCM。

对DM_BCM采用四阶龙格库塔法获得共振Cycler轨道SM_BCM，该SM_BCM轨道如图3所示。图中，各点反应了五个周期里航天器P在地月质心惯性坐标系O-XYZ下的位置。所述四阶龙格库塔法参考2006年发表的《龙格库塔法及其Mathematica实现》，第18卷，第2期，第72至第73页，作者陈誌敏。

在本发明中，采用多重打靶法对初始化的M_L模型进行修正得到修正后的同宿Cycler轨道动力学模型DM_L。所述DM_L是以图形(如图3A所示)进行表征的。

在一个周期的M_L模型上取500个等时间间隔的样点，任意一样点记为ini_d(i)(i＝1，…，500)；根据M_L模型采用四阶龙格库塔法积分，得到任意一点的位置速度状态量记为ini_c(i)(i＝1，…，500)；ini_c(i，：)是一个500×6的矩阵。

在每个时间间隔上，选取100个样点，并对样点采用ODE45积分，取出末点的位置速度状态量θ(ini_c(i，：))。记末点的位置速度状态量θ(ini_c(i，：))与任意一点的位置速度状态量ini_c(i，：)的差值为G(i，：)，即G(i，：)=θ(ini_c(i，：))-ini_c(i，：)。应用差值G(i，：)是使得解算值和实际计算值能够无误差，在仿真过程中应尽量将误差缩减到最小。

对于上述差值G(i，：)采用牛顿迭代方法求极小值，在设定精度为1×10^-10条件下，经过8次迭代，达到设定精度。

第1次迭代	1.4925020681966167×10-1
		第2次迭代	8.5043845337684230×10-2
第3次迭代	2.3461313323453772×10-2
		第4次迭代	1.9075592366893615×10-3

第5次迭代	5.5915036547800547×10-4
		第6次迭代	8.9114511750870658×10-7
第7次迭代	4.5782899616269825×10-9
		第8次迭代	5.1435802854146783×10-10

若大于所设定的精度，继续以牛顿迭代法重复上述迭代，直至小于所设定的精度为止。由此，可以得到修正后的模型DM_L。

对DM_L采用四阶龙格库塔法获得共振Cycler轨道SM_L，该SM_L轨道如图3A所示。图中，各点反应了航天器P在地月质心惯性坐标系O—XYZ下的位置。

对DM_L采用四阶龙格库塔法获得同宿Cycler轨道SM_L。

(四)Lambert转移轨道获取模块40

Lambert问题是航天动力学中的两点边界值问题，在航天器交会、导弹拦截、星际航行等领域中有广泛的应用。如图4所示，航天器P的初始端点E₁、终端E₂的位置矢量分别为L₁和L₂，椭圆转移轨道的焦点位于地心P₁。初始端点E₁的时间记为t₁，终端E₂的时间记为t₂，θ为转移角。

根据Lambert飞行时间定理，椭圆转移轨道上任意两点之间的转移时间与椭圆转移轨道的长半轴ra、两点半径之和(L₁+L₂)以及中心角θ有关，则表示为：

t_f=W(ra，(L₁+L₂)，R_P)   (4)

若L₁与L₂之和为常数，长半轴ra为常数，初始端点E₁与终端E₂之间的距离R_P为常数，则从初始端点E₁至终端E₂的飞行转移时间t_f也是常数。椭圆转移轨道的确定以及两点间速度的选择是Lambert问题的关键，它被描述为如下的高斯问题：追踪航天器P出发处位置矢量(L₁)与速度矢量(v₁)，航天器P终端位置矢量(L₂)与速度矢量(v₂)，飞行转移时间为t_f，航天器P在椭圆转移轨道E₁点处的初始速度为v₁₀，航天器P在椭圆转移轨道E₂点处的结束速度为v₂₀。该高斯问题的目标是为了求解初始位置速度增量Δv₁与终端位置速度增量Δv₂，进而求的施加的脉冲推力大小。

对于高斯问题可以由如下的超越方程组求得：

L₂=k×L₁+g×v₁₀   (5)

$v_{20}=\stackrel{\·}{k}\×L_1+\stackrel{\·}{g}\×v_{10}---\left(6\right)$

拉格朗日系数一 $k=1-\frac{{\mathrm{\μL}}_2}{h^2}(1-\cos \mathrm{\θ});$

拉格朗日系数二 $\stackrel{\·}{k}=\frac{\mathrm{\μ}}{h}\frac{1-\cos \mathrm{\θ}}{\sin \mathrm{\θ}}\[\frac{\mathrm{\μ}}{h}(1-\cos \mathrm{\θ})-\frac{1}{L_1}-\frac{1}{L_2}\];$

拉格朗日系数三

拉格朗日系数四

μ为次天体与主天体质量之比，在CR3BP模型中，取值为0.01215；h为航天器P在不同位置处的椭圆转移轨道上的变量。

在公式(5)、公式(6)中，已知L₁和v₁₀，或者L₂和v₂₀就可以确定出航天器P运行的椭圆转移轨道。显然，一旦确定了拉格朗日系数k、g、高斯问题便可迎刃而解。

在本发明中，对于公式(5)、公式(6)采用全局变量算法进行求解。

所述高斯-全局变量复合算法参考2008年出版的《空间交会对接任务规划》，第81页至第100页内容，第142页至第144页，作者唐国金等。

假设月球至共振Cycler轨道转移时间为3天，Lambert转移轨道起点距月球表面高度100km。共振Cycler轨道至地球转移时间为0.5天，Lambert转移轨道终点距地球表面高度100km。假设地球至共振Cycler轨道转移时间为0.5天，Lambert转移轨道起点距地球表面高度100km。共振Cycler轨道至月球转移时间为1天，Lambert转移轨道终点距月地球表面高度100km。

采用高斯-全局变量复合算法对SM_BCM进行处理，得到第一套Lambert转移轨道 $S_{\mathrm{BCM}}=\{s_1^{\mathrm{BCM}},s_2^{\mathrm{BCM}},s_3^{\mathrm{BCM}},s_4^{\mathrm{BCM}}\},$ 相应的Lambert轨道仿真结果如图5所示。图中，在地月质心惯性坐标系O-XYZ下有4段虚线，第一条虚线由Cycler轨道至月球停泊轨道的Lambert转移轨道；第二条虚线

由月球停泊轨道至Cycler轨道的Lambert转移轨道；第三条虚线

由地球停泊轨道至Cycler轨道的Lambert转移轨道；第四条虚线由Cycler轨道至地球停泊轨道的Lambert转移轨道。实线为SM_BCM轨道。

经典二体Lambert问题是典型的两点边值问题，采用微分修正的方法进行修正第一套Lambert转移轨道

修正过程为：

根据L₁、L₂和θ可计算二体Lambert转移轨道。Lambert变轨策略给出位置交会的转移轨道，控制量为初始端点E₁位置处的脉冲速度增量，其方向角记为β；微分修正算法将改进初始位置速度增量Δv₁及转移时间t_f，以实现终端点E₂位置处的交会。贴近真值的迭代初值，可以保证微分修正算法迭代过程的收敛。为了将航天器P导引到目标位置L₂(即理论终端位置矢量)，每次迭代过程都将轨道积分至L₂处，而t_f即取为该轨道积分时间。显然，初始端点E₁位置处的轨道速度v₁₀的变化Δv₁将导致轨道积分时间t_f的变化，记为Δt_f。考查第m次迭代(迭代m次的前一次记为m-1，迭代m次的后一次记为m+1，则迭代m次记为当前次)，对终端点E₂位置处的位置矢量作一阶Taylor展开，可得：

$L(t_f+\mathrm{\Δ}{t}_f,v_{10}+\mathrm{\Δ}{v}_1)=L_2+\frac{{\∂L}_2}{{\∂v}_{10}}\×\mathrm{\Δ}{v}_1+v_{20}\×\mathrm{\Δ}{t}_f---\left(7\right)$

L(t_f+Δt_f，v₁₀+Δv₁)表示航天器P终端位置实际矢量关于转移时间、在椭圆转移轨道E₁点处的初始速度之间的函数；

表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数，简化为 $M=\frac{{\∂L}_2}{{\∂v}_{10}};$

在本发明中，速度修正量Δv₁应使得：

L(t_f+Δt_f，v₁₀+Δv₁)＝L₂    (8)

整理式(7)和式(8)，可得：

$\mathrm{\δ}{L}_2=L_2-\stackrel{\‾}{L_2}=\frac{{\∂L}_2}{{\∂v}_{10}}\×\mathrm{\Δ}{v}_1-v_{20}\×\mathrm{\Δ}{t}_f---\left(9\right)$

表示航天器P的实际终端位置矢量。

不妨取

和L₂具有相同的横坐标分量，在日地旋转坐标系O_S-X_SY_SZ_S下，并记终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数

展开式(9)，可得：

$\left[\begin{array}{c}0\\ \mathrm{\δy}\\ \mathrm{\δz}\end{array}\right]=-\left[\begin{array}{cccc}{M}_{11}& M_{12}& M_{13}& v_{20}^x\\ M_{21}& M_{22}& M_{23}& v_{20}^y\\ M_{31}& M_{32}& M_{33}& v_{20}^z\end{array}\right]\×\left[\begin{array}{c}\mathrm{\Δ}{v}_1^x\\ \mathrm{\Δ}{v}_1^y\\ {\mathrm{\Δv}}_1^z\\ {\mathrm{\Δt}}_f\end{array}\right]---\left(10\right)$

δy表示日地旋转坐标系O_S-X_SY_SZ_S下的Y_S轴的位置差值；

δz表示日地旋转坐标系O_S-X_SY_SZ_S下的Z_S轴的位置差值；

M₁₁表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第一行第一列元素；

M₁₂表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第一行第二列元素；

M₁₃表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第一行第三列元素；

M₂₁表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第二行第一列元素；

M₂₂表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第二行第二列元素；

M₂₃表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第二行第三列元素；

M₃₁表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第三行第一列元素；

M₃₂表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第三行第二列元素；

M₃₃表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第三行第三列元素；

表示终端位置E₂处结束速度在X_S轴的速度分量；

表示终端位置E₂处结束速度在Y_S轴的速度分量；

表示终端位置E₂处结束速度在Z_S轴的速度分量；

表示初始端点E₁位置处初始速度增量在X_S轴的速度分量；

表示初始端点E₁位置处初始速度增量在Y_S轴的速度分量；

表示初始端点E₁位置处初始速度增量在Z_S轴的速度分量；

Δt_f表示初始端点E₁位置处的轨道速度v₁₀的增量Δv₁导致的轨道积分时间的变化。

令 $C=\left[\begin{array}{cccc}{M}_{11}& M_{12}& M_{13}& v_{20}^x\\ M_{21}& M_{22}& M_{23}& v_{20}^y\\ M_{31}& M_{32}& M_{33}& v_{20}^z\end{array}\right],$ 则有：

$\left[\begin{array}{c}{\mathrm{\Δv}}_1^x\\ {\mathrm{\Δv}}_1^y\\ {\mathrm{\Δv}}_1^z\\ {\mathrm{\Δt}}_f\end{array}\right]=-{\left(C^{T}C\right)}^{-1}{C}^T\left[\begin{array}{c}0\\ \mathrm{\δy}\\ \mathrm{\δz}\end{array}\right]---\left(11\right)$

C^T表示矩阵C的倒置。

故相对于前一次m-1的速度来讲，则当前次m的速度增量表征为：

$\mathrm{\Δ}{v}_1=\left[\begin{array}{ccc}{\mathrm{\Δv}}_1^x& {\mathrm{\Δv}}_1^y& {\mathrm{\Δv}}_1^z\end{array}\right]---\left(12\right)$

在设定精度后，反复迭代，直至满足精度。针对月球停泊轨道起飞至共振型Cycler轨道的转移轨道，设定精度为1×10^-8迭代结果如下表所示。

第1次迭代	2.0620622902499468×100
		第2次迭代	2.7641929155045797×10-1
第3次迭代	7.6981072705594636×10-2
		第4次迭代	8.3862124467837824×10-3
第5次迭代	1.7139690673584045×10-4
		第6次迭代	1.0442627289698014×10-7

经过6次迭代，达到设定精度。

针对共振型Cycler轨道至月球停泊轨道的转移轨道，设定精度为1×10^-8迭代结果如下表所示。

第1次迭代	2.5421929920011647×100
		第2次迭代	1.2546155043467997×10-1
第3次迭代	6.5637444323567223×10-2
		第4次迭代	5.3832424435867889×10-3
第5次迭代	2.2263784592810393×10-4
		第6次迭代	2.0143221442658014×10-7

经过6次迭代，达到设定精度。

针对地球停泊轨道起飞至共振型Cycler轨道的转移轨道，设定精度为1×10^-8迭代结果如下表所示。

第1次迭代	8.5792468327776916×10-1
		第2次迭代	7.9205098920078354×10-1
第3次迭代	1.4596634587845803×10-1
		第4次迭代	2.5378649526028702×10-3
第5次迭代	7.2114508959175548×10-7

经过5次迭代，达到设定精度。

针对共振型Cycler轨道至地球停泊轨道的转移轨道，设定精度为1×10^-8迭代结果如下表所示。

第1次迭代	8.2182773945012833×10-1
		第2次迭代	6.8909282035057489×10-1
第3次迭代	2.4658347859604583×10-1
		第4次迭代	3.0359578642687003×10-3

第5次迭代	1.2278451094586332×10-4
		第6次迭代	5.2845269283793103×10-5
第7次迭代	6.0899512115445758×10-7

经过7次迭代，达到设定精度。

在本发明中，采用微分修正算法对 $S_{\mathrm{BCM}}=\{s_1^{\mathrm{BCM}},s_2^{\mathrm{BCM}},s_3^{\mathrm{BCM}},s_4^{\mathrm{BCM}}\};$ 进行修正，得到第一套修正后的Lambert转移轨道 ${\mathrm{DS}}_{\mathrm{BCM}}=\{{\mathrm{ds}}_1^{\mathrm{BCM}},{\mathrm{ds}}_2^{\mathrm{BCM}},{\mathrm{ds}}_3^{\mathrm{BCM}},{\mathrm{ds}}_4^{\mathrm{BCM}}\},$ 如图5A、图5B所示：图中，修正后为粗实线。航天器在实际转移轨道上的飞行情况更接近粗实线，更符合运动规律。对比图5与图5A，其中

为

的修正，

为的修正。对比图5与图5B，其中

为

的修正，

为

的修正。

(一)分析航天器由地面起飞至Cycler轨道的交会

首先进行发射窗口分析：

白赤(地球赤道面)交角呈周月变化，2005年3月21日至2006年3月21日期间的变化范围在[0°，η]，其中η取值在28.4°～28.7°之间变化。西昌发射纬度为27.9°，文昌发射纬度为19°，则两者均存在发射机会可将空间站直接射入白道面，且西昌直接入轨的窗口小于文昌；如图6所示，横坐标为时间，纵坐标为白赤(地球赤道面)交角。对于文昌，发射窗口的宽度约为7天，间隔约为7天；对于西昌，发射窗口的宽度约为2天，间隔约为13天。在发射窗口内，应选择合适的发射时机(1天1次)以将升交点赤经捕获到白道面内。

西昌的发射纬度参考2007年发表的《“火龙出水”的西昌发射场——素描中国四大卫星发射中心(1)》，第17页，作者晓波。

文昌的发射纬度参考2012年发表的《海南文昌火箭发射场雷电环境分析》，第183页，作者高燚等。

然后进行交会窗口分析：

空间站由停泊轨道(高度为100km)转移到共振Cycler轨道，采取霍曼转移方式或Lambert转移方式进行燃料消耗量与转移时间的分析，如图6A和图6B所示。图6A中为最省燃料霍曼转移所需速度增量为3192m/s，转移时间为7天。图6B中为最费燃料霍曼转移所需速度增量为4093.5m/s，转移时间为66.5min。

由于图6A和图6B达到远地点的时间大致相等(约为7天)，而图6A燃料更为节省，将作为标称交会轨迹：记标称交会为第1天发射轨迹，若因为标称交会发射延误，选择第2、3、4、5、6、7天进行发射，则所需燃料依次由3192km／s增加至4093.5m／s；记标称交会为第7天发射轨迹，也可选择从倒数天数记，-6、-5、-4、-3、-2、-1天进行发射，则所需燃料依次由3192km／s增加至4093.5m／s；记标称交会为第1～7天发射轨迹，则所需燃料依次由3192km／s增加至4093.5m／s。连续7天的发射轨迹如图6C和图6D所示，在地月质心惯性坐标系O-XYZ下发射轨迹如图6C所示，地月质心旋转坐标系O-xyz下发射轨迹如图6D所示。

(二)分析航天器由月面起飞至共振Cycler轨道的交会

首先进行发射窗口分析：

白赤(月球赤道面)交角呈周月变化，2005年3月21日至2006年3月21日期间的变化范围在[0°，ε]，其中ε取值在6.64°～6.86°之间变化，如图7所示，横坐标为时间，纵坐标为白赤(月球赤道面)交角。变化周期约为13.6天。对于纬度高于6.64°的着陆点，月面起飞后需要调整轨道平面，方可进入共振Cycler轨道。以南北极着陆点为例，月面起飞进入100km环月极轨；共振Cycler轨道捕获过程如下：抬升远月点至20000km，所需速度增量约为577m／s；远月点进行倾角修正，所需速度增量约为287.3m／s。远月点的抬升可为下一步共振Cycler轨道交会节省部分能量。月球公转周期和自转周期相同，每个月只有两次机会可将升交点赤经捕获到白道面内。

然后进行交会窗口分析：

以交会时间和共振Cycler轨道相位为变量，以月面停泊轨道高度100km作为约束条件，搜索速度增量小于等于1500m／s的交会窗口。地理经度取决于升交点赤经，根据“停泊轨道高度100km”作为约束条件的搜索结果如图7A所示，在图7A中，T1点、T2点、T3点和T4点表示符合搜索条件、且最优的交会情况。如图7B所示的地理经度可分布在[-170°，-135°]∪[-43°，-5°]∪[45°，108°]∪[150°，165°]。在上述区间内的着陆点，通过月面起飞后有望直接进入白道面内；上述区间外的着陆点需要等待半个月进行升交点赤经捕获。因此，若选择指定区域内的着陆点，月面起飞后即进入白道面内：对于共振Cycler轨道相位[0.05,0.25]∪[0.8,1]内可随时起飞交会。

在本发明中，采用高斯-全局变量复合算法对SM_L进行处理，得到第二套Lambert转移轨道

对SM_L的处理方法与对SM_BCM的处理方法是相同的。为了区分说明，在引述的关系式上添加字母H进行区别。

根据Lambert飞行时间定理，椭圆转移轨道上任意两点之间的转移时间与椭圆转移轨道的长半轴Hra、两点半径之和(HL₁+HL₂)以及中心角Hθ有关，则表示为：

Ht_f=W(Hra,(HL₁+HL₂)，HR_P)             (13)

若HL₁与HL₂之和为常数，长半轴Hra为常数，初始端点E₁与终端E₂之间的距离HR_P为常数，则从初始端点E₁至终端E₂的飞行转移时间Ht_f也是常数。椭圆转移轨道的确定以及两点间速度的选择是Lambert问题的关键，它被描述为如下的高斯问题：追踪航天器P出发处位置矢量(HL₁)与速度矢量(HV₁)，航天器P终端位置矢量(HL₂)与速度矢量(Hv₂)，飞行转移时间为Ht_f，航天器P在椭圆转移轨道E₁点处的初始速度为Hv₁₀，航天器P在椭圆转移轨道E₂点处的结束速度为Hv₂₀。该高斯问题的目标是为了求解初始位置速度增量△Hv₁与终端位置速度增量ΔHv₂，进而求的施加的脉冲推力大小。

在本发明中，高斯问题可以由如下的超越方程组求得：

HL₂=k×HL₁+g×Hv₁₀      (14)

${\mathrm{Hv}}_{20}=\stackrel{\·}{k}\×{\mathrm{HL}}_1+\stackrel{\·}{g}\×{\mathrm{Hv}}_{10}---\left(15\right)$

拉格朗日系数一 $k=1-\frac{\mathrm{\μ}\×{\mathrm{HL}}_2}{{\mathrm{Hh}}^2}(1-\cos \mathrm{H\θ});$

拉格朗日系数二 $\stackrel{\·}{k}=\frac{\mathrm{\μ}}{\mathrm{Hh}}\frac{1-\cos \mathrm{H\θ}}{\sin \mathrm{H\θ}}[\frac{\mathrm{\μ}}{\mathrm{Hh}}(1-\cos \mathrm{H\θ})-\frac{1}{{\mathrm{HL}}_1}-\frac{1}{{\mathrm{HL}}_2}];$

拉格朗日系数三 $g=\frac{{\mathrm{HL}}_1\×{\mathrm{HL}}_2\×\sin \mathrm{H\θ}}{\mathrm{Hh}};$

拉格朗日系数四

μ为次天体与主天体质量之比，在BCM模型中，取值为0.000003003；h为航天器P在不同位置处的椭圆转移轨道上的变量。

在公式(14)、公式(15)中，已知HL₁和Hv₁₀，或者HL₂和Hv₂₀就可以确定出航天器P运行的椭圆转移轨道。显然，一旦确定了拉格朗日系数k、g、高斯问题便可迎刃而解。

在本发明中，对于公式(14)、公式(15)采用全局变量算法进行求解。

假设地球至同宿轨道转移时间为1天，转移轨道起点距地球表面高度100km。同宿轨道至月球转移时间为4天，转移轨道终点距月球表面高度100km。根据二体Lambert问题解法，相应的同宿Lambert轨道仿真结果如图7C所示。

采用微分修正算法对

进行修正，得到第二套修正后的Lambert转移轨道

第二套修正后转移轨道仿真结果如7D所示。

根据HL₁、HL₂和Hθ可计算二体Lambert转移轨道。Lambert变轨策略给出位置交会的转移轨道，控制量为初始端点E₁位置处的脉冲速度增量，其方向角记为Hβ；微分修正算法将改进初始位置速度增量△Hv₁及转移时间Ht_f，以实现终端点E₂位置处的交会。贴近真值的迭代初值，可以保证微分修正算法迭代过程的收敛。为了将航天器P导引到目标位置L₂(即理论终端位置矢量)，每次迭代过程都将轨道积分至HL₂处，而Ht_f即取为该轨道积分时间。显然，初始端点E₁位置处的轨道速度Hv₁₀的变化ΔHv₁将导致轨道积分时间Ht_f的变化，记为ΔHt_f。考查第m次迭代(迭代m次的前一次记为m-1，迭代m次的后一次记为m+1，则迭代m次记为当前次)，对终端点E₂位置处的位置矢量作一阶Taylor展开，可得：

$\mathrm{HL}({\mathrm{Ht}}_f+\mathrm{\Δ}{\mathrm{Ht}}_f,{\mathrm{Hv}}_{10}+\mathrm{\Δ}{\mathrm{Hv}}_1)={\mathrm{HL}}_2+\frac{\∂{\mathrm{HL}}_2}{\∂{\mathrm{Hv}}_{10}}\×\mathrm{\Δ}{\mathrm{Hv}}_1+{\mathrm{Hv}}_{20}\×\mathrm{\Δ}{\mathrm{Ht}}_f---\left(16\right)$

HL(Ht_f+ΔHt_f，Hv₁₀+ΔHv₁)表示航天器P终端位置实际矢量关于转移时间、在椭圆转移轨道E₁点处的初始速度之间的函数；

表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数，简化为 $\mathrm{HM}=\frac{\∂{\mathrm{HL}}_2}{\∂{\mathrm{Hv}}_{10}};$

在本发明中，速度修正量ΔHv₁应使得：

HL(Ht_f+ΔHt_f，Hv₁₀+ΔHv₁)=HL₂              (17)

整理式(16)和式(17)，可得：

$\mathrm{\δ}{\mathrm{HL}}_2={\mathrm{HL}}_2-{\stackrel{\‾}{\mathrm{HL}}}_2=\frac{\∂{\mathrm{HL}}_2}{\∂{\mathrm{Hv}}_{10}}\×\mathrm{\Δ}{\mathrm{Hv}}_1-{\mathrm{Hv}}_{20}\×\mathrm{\Δ}{\mathrm{Ht}}_f---\left(18\right)$

表示航天器P的实际终端位置矢量。

不妨取

和L₂具有相同的横坐标分量，在日地旋转坐标系O_S-X_SY_SZ_S下，并记终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数

展开式(18)，可得：

$\left[\begin{array}{c}0\\ \mathrm{\δHy}\\ \mathrm{\δHz}\end{array}\right]=-\left[\begin{array}{cccc}{\mathrm{HM}}_{11}& {\mathrm{HM}}_{12}& {\mathrm{HM}}_{13}& H{v}_{20}^x\\ {\mathrm{HM}}_{21}& {\mathrm{HM}}_{22}& {\mathrm{HM}}_{23}& H{v}_{20}^y\\ {\mathrm{HM}}_{31}& {\mathrm{HM}}_{32}& {\mathrm{HM}}_{33}& H{v}_{20}^z\end{array}\right]\×\left[\begin{array}{c}\mathrm{\ΔH}{v}_1^x\\ \mathrm{\ΔH}{v}_1^y\\ \mathrm{\ΔH}{v}_1^z\\ \mathrm{\ΔH}{t}_f\end{array}\right]---\left(19\right)$

δHy表示日地旋转坐标系O_S-X_SY_SZ_S下的Y_S轴的位置差值；

δHz表示日地旋转坐标系O_S-X_SY_SZ_S下的Z_S轴的位置差值；

HM₁₁表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第一行第一列元素；

HM₁₂表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第一行第二列元素；

HM₁₃表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第一行第三列元素；

HM₂₁表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第二行第一列元素；

HM₂₂表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第二行第二列元素；

HM₂₃表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第二行第三列元素；

HM₃₁表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第三行第一列元素；

HM₃₂表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第三行第二列元素；

HM₃₃表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第三行第三列元素；

表示终端位置E₂处结束速度在X_S轴的速度分量；

表示终端位置E₂处结束速度在Y_S轴的速度分量；

表示终端位置E₂处结束速度在Z_S轴的速度分量；

表示初始端点E₁位置处初始速度增量在X_S轴的速度分量；

表示初始端点E₁位置处初始速度增量在Y_S轴的速度分量；

表示初始端点E₁位置处初始速度增量在Z_S轴的速度分量；

△Ht_f表示初始端点E₁位置处的轨道速度v₁₀的增量△v₁导致的轨道积分时间的变化。

令 $\mathrm{HC}=\left[\begin{array}{cccc}{\mathrm{HM}}_{11}& {\mathrm{HM}}_{12}& {\mathrm{HM}}_{13}& H{v}_{20}^x\\ {\mathrm{HM}}_{21}& {\mathrm{HM}}_{22}& {\mathrm{HM}}_{23}& H{v}_{20}^y\\ {\mathrm{HM}}_{31}& {\mathrm{HM}}_{32}& {\mathrm{HM}}_{33}& H{v}_{20}^z\end{array}\right],$ 则有：

$\left[\begin{array}{c}\mathrm{\ΔH}{v}_1^x\\ \mathrm{\ΔH}{v}_1^y\\ \mathrm{\ΔH}{v}_1^z\\ \mathrm{\ΔH}{t}_f\end{array}\right]=-{\left({\mathrm{HC}}^T\mathrm{HC}\right)}^{-1}{\mathrm{HC}}^T\left[\begin{array}{c}0\\ \mathrm{\δHy}\\ \mathrm{\δHz}\end{array}\right]---\left(20\right)$

HC^T表示矩阵HC的倒置。

故相对于前一次m-1的速度来讲，则当前次m的速度增量表征为：

$\mathrm{\ΔH}{v}_1=\left[\begin{array}{ccc}\mathrm{\ΔH}{v}_1^x& \mathrm{\ΔH}{v}_1^y& \mathrm{\ΔH}{v}_1^z\end{array}\right]---\left(21\right)$

在设定精度后，反复迭代，直至满足精度。针对月球停泊轨道起飞至共振型Cycler轨道的转移轨道，设定精度为1×10^-15迭代结果如下表所示。

第1次迭代	2.8639968709823611×100
		第2次迭代	1.2827923060663518×100
第3次迭代	5.6581823976188161×10-1
		第4次迭代	2.1973874406340671×10-1
第5次迭代	1.5468617860521472×10-1
		第6次迭代	4.9764011654783211×10-2
第7次迭代	5.2642496344864611×10-3
		第8次迭代	2.0846898885215647×10-5
第9次迭代	2.2579833858870749×10-10
		第10次迭代	6.5486562934038266×10-14

经过10次迭代，达到设定精度。

针对地球停泊轨道至同宿型Cycler轨道的转移轨道，设定精度为1×10^-11迭代结果如下表所示。

第1次迭代	4.0437979717451373×100
		第2次迭代	2.1369334069815822×100
第3次迭代	9.8469011235887069×10-2
		第4次迭代	5.1247925851644960×10-4
第5次迭代	5.5164314416937229×10-8
		第6次迭代	3.3282881768155430×10-12

经过6次迭代，达到设定精度。

在本发明中，比较基于共振型Cycler轨道SM_BCM及同宿型Cycler轨道SM_L的地月转移任务的飞行时间(包括变轨时间)、燃料消耗。以每个任务的两次Lambert变轨的总速度冲量表征燃料消耗量。

1)共振型Cycler轨道方案：

由地球停泊轨道至月球停泊轨道的总飞行时间为6.014天。

在此期间的两次Lambert变轨的总速度冲量为4.626km/s。

2)同宿型Cycler轨道方案：

由地球停泊轨道至月球停泊轨道的总飞行时间为20.695天。

在此期间的两次Lambert变轨的总速度冲量为4.450km/s。

通过对比总飞行时间、总燃料消耗，可以看出，基于Cycler轨道的地月转移方案，总飞行时间较长，总燃料消耗较少。这说明，基于Cycler轨道的地月转移方案的设想是可行的。两类Cycler轨道具有良好的周期性。共振型Cycler轨道的周期可以设定，航天器可以根据任务需要，有计划地往返于地月之间，规律的进行物资、人员的传输。同宿型Cycler轨道的周期也是可以确定的，并且经过了地月系LL₁点，在完成地月往返运输任务的同时，还可以对地月系LL₁点处布置的科学设备进行更换和检修。基于这两类Cycler轨道的地月转移方案飞行时间较长，总燃料消耗较少。合理利用其周期性好的优势，可以设计出地月往返“公交”***，对于小推力航天器的深空探测，地月***物资传输有重大的意义。

Claims (8)

1.一种基于Cycler轨道的地月往返的任务仿真***，其特征在于：该仿真***包括有构建共振Cycler轨道模型(10)、构建同宿Cycler轨道模型(20)、多重打靶法轨道修正模块(30)和Lambert转移轨道获取模块(40)；

所述构建共振Cycler轨道模型(10)第一方面依据CR3BP模型构建M_CR3BP动力学模型；第二方面采用BCM模型对所述的M_CR3BP动力学模型进行优化，得到M_BCM动力学模型；

所述构建同宿Cycler轨道模型(20)依据BCM模型构建M_L动力学模型；

所述多重打靶法轨道修正模块(30)第一方面采用多重打靶法分别对M_BCM进行修正，得到修正后的共振Cycler轨道动力学模型DM_BCM；第二方面对DM_BCM采用四阶龙格库塔法获得共振Cycler轨道SM_BCM；第三方面采用多重打靶法分别对M_L进行修正，得到修正后的同宿Cycler轨道动力学模型DM_L；第四方面对DM_L采用四阶龙格库塔法获得同宿Cycler轨道SM_L；

所述Lambert转移轨道获取模块(40)第一方面采用高斯-全局变量复合算法对SM_BCM进行处理，得到第一套Lambert转移轨道 $S_{\mathrm{BCM}}=\{s_1^{\mathrm{BCM}},s_2^{\mathrm{BCM}},s_3^{\mathrm{BCM}},s_4^{\mathrm{BCM}}\};$ 第二方面采用微分修正算法对

进行修正，得到第一套修正后的Lambert转移轨道 ${\mathrm{DS}}_{\mathrm{BCM}}=\{{\mathrm{ds}}_1^{\mathrm{BCM}},{\mathrm{ds}}_2^{\mathrm{BCM}},{\mathrm{ds}}_3^{\mathrm{BCM}},{\mathrm{ds}}_4^{\mathrm{BCM}}\};$ 第三方面依据发射点纬度和白赤交角变化规律对SM_BCM进行发射窗口及交会窗口的分析，获取航天器发射和入轨的时机；第四方面采用高斯-全局变量复合算法对SM_L进行处理，得到第二套Lambert转移轨道

第五方面采用微分修正算法对

进行修正，得到第二套修正后的Lambert转移轨道

第六方面通过比较地月往返任务的燃料消耗量、总飞行时间能够获知SM_BCM和SM_L轨道各自的优势，为小推力航天器的深空探测的地月往返任务提供优化设计指标。

2.根据权利要求1所述的基于Cycler轨道的地月往返的任务仿真***，其特征在于：所述多重打靶法轨道修正模块(30)中，在一个周期的M_BCM模型上取500个等时间间隔的样点，任意一样点记为ini_b(i)(i=1，…，500)；根据M_BCM模型采用四阶龙格库塔法积分，得到任意一点的位置速度状态量记为in_a(i)(i=1，…，500)；ini_a(i，：)是一个500×6的矩阵；在每个时间间隔上，选取100个样点，并对样点采用ODE45积分，取出末点的位置速度状态量φ(ini_a(i，：))；记末点的位置速度状态量φ(ini_a(i，：))与任意一点的位置速度状态量ini_a(i，：)的差值为F(i，：)，即F(i，：)=φ(ini_a(i，：))-ini_a(i，：)；应用差值F(i，：)是使得解算值和实际计算值能够无误差，在仿真过程中应尽量将误差缩减到最小；对于上述差值F(i，：)采用牛顿迭代方法求极小值。

3.根据权利要求1所述的基于Cycler轨道的地月往返的任务仿真***，其特征在于：所述Lambert转移轨道获取模块(40)中，根据L₁、L₂和θ可计算二体Lambert转移轨道；Lambert变轨策略给出位置交会的转移轨道，控制量为初始端点E₁位置处的脉冲速度增量，其方向角记为β；微分修正算法将改进初始位置速度增量Δv₁及转移时间t_f，以实现终端点E₂位置处的交会；为了将航天器P导引到目标位置L₂，每次迭代过程都将轨道积分至L₂处，而t_f即取为该轨道积分时间；显然，初始端点E₁位置处的轨道速度v₁₀的变化Δv₁将导致轨道积分时间t_f的变化，记为Δt_f；考查第m次迭代，对终端点E₂位置处的位置矢量作一阶Taylor展开，可得 $L(t_f+\mathrm{\Δ}{t}_f,v_{10}+{\mathrm{\Δv}}_1)=L_2+\frac{{\∂L}_2}{{\∂v}_{10}}\×{\mathrm{\Δv}}_1+v_{20}\×{\mathrm{\Δt}}_f;$ L(t_f+Δt_f，v₁₀+Δv₁)表示航天器P终端位置实际矢量关于转移时间、在椭圆转移轨道E₁点处的初始速度之间的函数；

表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数，简化为

速度修正量Δv₁应使得L(t_f+Δt_f，v₁₀+Δv₁)=L₂成立，则有 ${\mathrm{\δL}}_2=L_2-\stackrel{\‾}{L_2}=\frac{{\∂L}_2}{{\∂v}_{10}}\×\mathrm{\Δ}{v}_1-v_{20}\×\mathrm{\Δ}{t}_f;$

表示航天器P的实际终端位置矢量；

取

和L₂具有相同的横坐标分量，在日地旋转坐标系O_S-X_SY_SZ_S下，并记终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数 $M=\frac{{\∂L}_2}{{\∂v}_{10}},$ 则有 $\left[\begin{array}{c}0\\ \mathrm{\δy}\\ \mathrm{\δz}\end{array}\right]=-\left[\begin{array}{cccc}{M}_{11}& M_{12}& M_{13}& v_{20}^x\\ M_{21}& M_{22}& M_{23}& v_{20}^y\\ M_{31}& M_{32}& M_{33}& v_{20}^z\end{array}\right]\×\left[\begin{array}{c}{\mathrm{\Δv}}_1^x\\ {\mathrm{\Δv}}_1^y\\ {\mathrm{\Δv}}_1^z\\ {\mathrm{\Δt}}_f\end{array}\right];$

δy表示日地旋转坐标系O_S-X_SY_SZ_S下的Y_S轴的位置差值；

δz表示日地旋转坐标系O_S-X_SY_SZ_S下的Z_S轴的位置差值；

M₁₁表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第一行第一列元素；

M₁₂表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第一行第二列元素；

M₁₃表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第一行第三列元素；

M₂₁表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第二行第一列元素；

M₂₂表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第二行第二列元素；

M₂₃表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第二行第三列元素；

M₃₁表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第三行第一列元素；

M₃₂表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第三行第二列元素；

M₃₃表示终端位置矢量与椭圆转移轨道E₁点处的初始速度之间的偏导数矩阵的第三行第三列元素；

表示终端位置E₂处结束速度在X_S轴的速度分量；

表示终端位置E₂处结束速度在Y_S轴的速度分量；

表示终端位置E₂处结束速度在Z_S轴的速度分量；

表示初始端点E₁位置处初始速度增量在X_S轴的速度分量；

表示初始端点E₁位置处初始速度增量在Y_S轴的速度分量；

表示初始端点E₁位置处初始速度增量在Z_S轴的速度分量；

Δt_f表示初始端点E₁位置处的轨道速度v₁₀的增量Δv₁导致的轨道积分时间的变化；

令 $C=\left[\begin{array}{cccc}{M}_{11}& M_{12}& M_{13}& v_{20}^x\\ M_{21}& M_{22}& M_{23}& v_{20}^y\\ M_{31}& M_{32}& M_{33}& v_{20}^z\end{array}\right],$ 则有 $\left[\begin{array}{c}{\mathrm{\Δv}}_1^x\\ {\mathrm{\Δv}}_1^y\\ {\mathrm{\Δv}}_1^z\\ {\mathrm{\Δt}}_f\end{array}\right]=-{\left(C^{T}C\right)}^{-1}{C}^T\left[\begin{array}{c}0\\ \mathrm{\δy}\\ \mathrm{\δz}\end{array}\right],$ C^T表示矩阵C的倒置；故相对于前一次m-1的速度来讲，则当前次m的速度增量表征为 ${\mathrm{\Δv}}_1=\left[\begin{array}{ccc}{\mathrm{\Δv}}_1^x& {\mathrm{\Δv}}_1^y& {\mathrm{\Δv}}_1^z\end{array}\right].$

4.根据权利要求1所述的基于Cycler轨道的地月往返的任务仿真***，其特征在于：依据CR3BP模型得到共振Cycler轨道动力学模型 $M_{\mathrm{CR}3\mathrm{BP}}=\left\{\begin{array}{c}\frac{\∂U}{\∂x}={\mathrm{HA}}^x-{2\mathrm{VA}}^y\\ \frac{\∂U}{\∂y}={\mathrm{HA}}^y+2{\mathrm{VA}}^x\\ \frac{\∂U}{\∂z}={\mathrm{HA}}^z\end{array}\right..$

5.根据权利要求1所述的基于Cycler轨道的地月往返的任务仿真***，其特征在于：依据BCM模型得到同宿Cycler轨道动力学模型 $M_L=\left\{\begin{array}{c}{\mathrm{HC}}^{X_S}={2\mathrm{VC}}^{Y_S}+{\mathrm{xs}}_P-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{xs}}_P+{\mathrm{\μ}}_2}{r_{\mathrm{PS}}^2}-\\ {\mathrm{\μ}}_2\frac{{\mathrm{xs}}_P-1+{\mathrm{\μ}}_2}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{xs}}_P-x_{P_2}}{r_{\mathrm{PM}}^3}\\ {\mathrm{HC}}^{Y_S}=-{2\mathrm{VC}}^{X_S}+{\mathrm{ys}}_P-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{ys}}_P}{r_{\mathrm{PS}}^3}-\\ {\mathrm{\μ}}_2\frac{{\mathrm{ys}}_P}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{ys}}_P-y_{P_2}}{r_{\mathrm{PM}}^3}\\ {\mathrm{HC}}^{Z_S}=-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{zs}}_P}{r_{\mathrm{PS}}^3}-{\mathrm{\μ}}_2\frac{{\mathrm{zs}}_P}{r_{\mathrm{PE}}^3}-m\frac{{\mathrm{zs}}_P}{r_{\mathrm{PM}}^3}\end{array}\right..$

6.根据权利要求1所述的基于Cycler轨道的地月往返的任务仿真***，其特征在于：利用发射点纬度和白赤交角变化规律所得，空间站由停泊轨道转移到共振Cycler轨道，采取霍曼转移方式或Lambert转移方式进行燃料消耗量与转移时间的分析，最省燃料霍曼转移所需速度增量为3192m／s，转移时间为7天；最费燃料霍曼转移所需速度增量为4093.5m／s，转移时间为66.5min；将最省燃料霍曼转移作为标称交会轨迹，若因为标称交会发射延误，选择第2、3、4、5、6、7天进行发射，则所需燃料依次由3192km／s增加至4093.5m／s；记标称交会为第7天发射轨迹，也可选择从倒数天数记，-6、-5、-4、-3、-2、-1天进行发射，则所需燃料依次由3192km／s增加至4093.5m／s。

7.根据权利要求1所述的基于Cycler轨道的地月往返的任务仿真***，其特征在于：利用发射点纬度和白赤交角变化规律所得，以交会时间和共振Cycler轨道相位为变量，以月面停泊轨道高度100km作为约束条件，搜索速度增量小于等于1500m／s的交会窗口；地理经度取决于升交点赤经，根据停泊轨道高度l00km作为约束条件能够搜索得到最优的交会情况；在上述区间内的着陆点，通过月面起飞后有望直接进入白道面内；上述区间外的着陆点需要等待半个月进行升交点赤经捕获。

8.根据权利要求1所述的基于Cycler轨道的地月往返的任务仿真***，其特征在于：比较基于共振型Cycler轨道SM_BCM及同宿型Cycler轨道SM_L的地月转移任务的飞行时间、燃料消耗；以每个任务的两次Lambert变轨的总速度冲量表征燃料消耗量；

所述共振型Cycler轨道：由地球停泊轨道至月球停泊轨道的总飞行时间为6.014天；在此期间的两次Lambert变轨的总速度冲量为4.626km／s；

所述同宿型Cycler轨道：由地球停泊轨道至月球停泊轨道的总飞行时间为20.695天；在此期间的两次Lambert变轨的总速度冲量为4.450km/s。

Images