一种在全太阳系星历模型中构建平动点轨道的数值方法

1.本发明涉及一种在全太阳系星历模型中构建平动点轨道的数值方法，属于轨道动力学技术领域


背景技术：

2.地月探测系统是否可以被建模为圆形限制性三体问题模型(简称crtbp)依赖于所关注的问题，在地月系统的影响范围外，真正的地月探测系统可以建模为crtbp模型，而在地月系统的影响范围内，月球的偏心率和其他行星对探测器的直接引力影响以及对地月系统的间接引力影响均不容忽视。因而地月平动点任务轨道设计需要一个高保真的模型和对外部扰动的更好理解。并且，在crtbp中的周期轨道并不存在于地月共线平动点附近。相反，考虑到两个中心流形和强的外部扰动，halo/lissayjous扰动轨道可能占主导地位。目前对于在整个太阳系引力模型中构造halo/lissayjous轨道的研究中，广泛用多步打靶法对其轨道进行求解，月球的偏心率和来自太阳系引力的影响会使halo/lissayjous轨道失去周期性，地球和平动点之间的距离不再是一个恒定的轨迹，而是沿着连接地球和月球的直线移动并且随时间变化的轨迹。由于时间跨度大，多步打靶的收敛性很大程度上取决于拼接点初始条件的准确性，因而获得多步打靶法中更为精确的拼接点成为在整个太阳系模型中构建halo/lissayjous轨道的主要问题。本专利技术所提出的一种在全太阳系模型中构建halo/lissayjous轨道改进的数值方法是基于poincar
é
截面的动力学行为分析，得到新的拼接点，利用改进的多步打靶法和序列二次规划构建halo/lissayjous轨道。


技术实现要素：

3.本发明的目的是建立地心会合坐标系下的动力学模型简化计算过程，并且获取更为准确的拼接点进行多步打靶，提出了一种在全太阳系星历模型中构建平动点轨道的数值方法。
4.本发明的目的是通过下述技术方案实现的：
5.所述的一种在全太阳系星历模型中构建平动点轨道的数值方法包含以下步骤：
6.步骤1:根据j2000下地月平动点附近探测器的运动方程推导出地心会合坐标系(grc)下地月平动点附近的轨道动力学方程
7.步骤2:根据上述推导的运动方程，定义两种类型的poincar
é
截面，利用其分析halo/lissayjous轨道动力学行为
8.步骤3:通过poincar
é
截面分析确定平动点瞬时会合坐标系(lisc)下crtbp模型中每段轨迹弧的初始条件
9.步骤4:确定平动点的瞬时状态
10.步骤5:基于平动点的瞬时状态确定地心会合坐标系(grc)下每段轨迹弧的初始条件,连接拼接点
11.步骤6:采用序列二次规划法得到所有内部拼接点位置速度的解并进行数值模拟
12.本发明所述的一种在全太阳系星历模型中构建平动点轨道的数值方法，对比已有技术，存在如下有益效果：
13.1.多步打靶法总是在惯性坐标系中使用一个星历模型来实现，在地心会合坐标系(grc)中提出了完整太阳系引力模型，并清楚地显示了grc相对于惯性坐标系的角速度，在这此模型中，所有天体轨道参数可以直接用星历计算出来，计算过程简便；
14.2.针对月球的偏心率和太阳系引力扰动导致探测器的不利偏移等实际因素，开发出一个高保真动力学模型,实用性强。
15.3.方法精度高，基于poincar
é
截面动力学行为分析搜索得到的多步打靶法中新的拼接点能够构造出更为精确的halo/lissayjous轨道，且方法适用性高，具有推广价值。
16.4.本专利技术能够处理考虑探测器具体几何尺寸和姿态信息而加入太阳辐射时的实际轨道设计任务，具有良好的工程价值。
附图说明
17.图1为本发明的在全太阳系星历模型中构建平动点轨道的数值方法流程图；
18.图2为本发明的坐标系统示意图；
19.图3为本发明的crtbp中lisc下halo轨道和lissayjous轨道的典型poincar
é
截面示意图；(a)为lisc下halo轨道的poincar
é
截面，(b)为lisc下lissayjous轨道的poincar
é
截面
20.图4为本发明的grc与l1sc下halo轨道的poincar
é
截面示意图；(a)为grc下halo轨道的poincar
é
截面，(b)为l1sc下halo轨道的poincar
é
截面
21.图5为本发明的grc与l1sc下lissayjous轨道的poincar
é
截面示意图；(a)为grc下lissayjous轨道的poincar
é
截面，(b)为l1sc下lissayjous轨道的poincar
é
截面
22.图6为本发明的grc与l2sc下halo轨道的poincar
é
截面示意图；(a)为grc下halo轨道的poincar
é
截面，(b)为l2sc下halo轨道的poincar
é
截面
23.图7为本发明的grc与l2sc下lissayjous轨道的poincar
é
截面示意图；(a)为grc下lissayjous轨道的poincar
é
截面，(b)为l2sc下lissayjous轨道的poincar
é
截面
24.图8为本发明的grc与l3sc下halo轨道的poincar
é
截面示意图；(a)为grc下halo轨道的poincar
é
截面，(b)为l3sc下halo轨道的poincar
é
截面
25.图9为本发明的在x-y平面上投影的地月系统的halo/lissayjous轨道上的不同相位角η值的位置示意图
26.图10为本发明的lisc和grc的几何关系示意图(a)rk与r
l1
的几何关系(b)rk与r
l2
的几何关系(c)rk与r
l3
的几何关系
27.图11(a)为本发明的crtbp下间断的初始猜测图
28.图11(b)为本发明的间断弧收敛到位置和速度处连续的解图
29.图12为本发明的crtbp模型中l1sc下的halo/lissayjous轨道示意图(a)为halo轨道，(b)为lissayjous轨道
30.图13为本发明的crtbp模型中l2sc下的halo/lissayjous轨道示意图(a)为halo轨道，(b)为lissayjous轨道
31.图14为本发明的crtbp模型中l3sc下的halo/lissayjous轨道示意图(a)为halo轨
道，(b)为lissayjous轨道
32.图15为本发明的拼接完成的l1点的halo扰动轨道示意图(a)grc下的轨迹(b)l1sc下的轨迹
33.图16为本发明的拼接完成的l1点的lissayjous扰动轨道示意图(a)grc下的轨迹(b)l1sc下的轨迹
34.图17为本发明的拼接完成的l2点的halo扰动轨道示意图(a)grc下的轨迹(b)l2sc下的轨迹
35.图18为本发明的拼接完成的l2点的lissayjous扰动轨道示意图(a)grc下的轨迹(b)l2sc下的轨迹
36.图19为本发明的拼接完成的l3点的halo扰动轨道示意图(a)grc下的轨迹(b)l3sc下的轨迹
具体实施方式
37.本发明提出了一种在全太阳系星历模型中构建平动点轨道的数值方法，主要包括以下步骤：
38.步骤1:根据j2000下地月平动点附近探测器的运动方程推导出地心会合坐标系(grc)下地月平动点附近的轨道动力学方程
39.图2所示，将平动点附近的运动定义在会合架构中。在整个太阳系引力模型中，航天器被看成是“无质量”的。
40.地心会合坐标系o-xyz(geocentric rotating coordinate system)，简称grc。原点在地心o，ox轴从地心指向月心，z轴垂直于地月系统的轨道平面并指向月球的瞬时轨道角动量方向，y轴与x轴、z轴构成右手直角坐标系。
41.地心赤道惯性坐标系o-xyz(j2000 geocentric equatorial coordinate system)，简称j2000。以地心o为原点，j2000.0历元时刻的地球平赤道面为参考面，x轴指向该历元时刻的平春分点，z轴与参考面的正法向一致，y轴与x轴、z轴构成右手直角坐标系。
42.为研究平动点周围的运动，将坐标系的原点移动到平动点；即li(i＝1,2,3)瞬时会合坐标系(lisynodic coordinate system)，简称lisc。l
i-xi、l
i-y
ii
、l
i-zi分别与ox轴、oy轴和oz轴平行，此坐标系是一个随时间变化的坐标系。
43.j2000下地月平动点附近探测器的运动方程为
[0044][0045]
其中μe，μm和μ
si
分别是地球、月球和太阳lisc和其他行星的引力常数，包括水星、金星、火星、木星、土星、天王星和海王星。r
p
是从地心到探测器的位置向量，rm和r
si
分别表示月球、太阳和其他行星相对于地心的位置向量。rm、r
p
、r
si
表示rm、r
p
和r
si
的标量形式
[0046][0047]
根据惯性坐标系与会合坐标系之间运动矢量的相对微分关系，得到以下关系：
[0048][0049]
其中，ω表示grc相对于j2000的角速度，和分别表示grc中从地球质心到的探测器的位置、速度和加速度
[0050]
由式(1)和式(3)，推导出grc下地月平动点附近探测器的运动方程
[0051][0052]
其中
[0053][0054][0055][0056]
根据grc的定义，获知到ω是月球相对于地球运动的轨道角速度，为得到ω的表达式，首先分析月球附近探测器的运动方程
[0057][0058]
其中r
msi
＝||r
m-r
si
||
[0059]
根据grc的定义，x轴、y轴、z轴方向的单位向量i、j、k可以表示为
[0060][0061]
其中月球角动量为
[0062][0063]
此外，单位向量i、j和k在空间上不是固定的，而是不断变化的方向；因此，它们的时间导数不是零。它们显然有一个恒定的大小(单位)，并且被附加在xyz框架上，它们都有相同的角速度ω。因此
[0064][0065]
因而角速度ω可以在grc表示为
[0066]
ω＝(ω
·
i)i+(ω
·
j)j+(ω
·
k)k
ꢀꢀꢀ
(12)
[0067]
根据式(11)和式(12)，可以得以下关系
[0068][0069]
将式(9)代入到式(13)中，角速度ω表达式变为
[0070][0071][0072]
其中
[0073][0074][0075]
将角速度ω表达式代入到式(4)中得到：
[0076][0077]
步骤2:根据上述推导的运动方程，定义两种类型的poincar
é
截面，利用其分析halo/lissayjous轨道动力学行为
[0078]
选择相同的历元时刻来映射poincar
é
截面，定义两种类型的poincar
é
截面，grc中
的x-y平面被定义为第一个poincar
é
截面来分析grc下halo/lissayjous轨道动力学行为；lisc中的x
i-yi平面被定义为第二个poincar
é
截面。
[0079]
l3点附近的周期轨道线性化平面内频率非常接近线性化的平面外频率，从而总呈现出halo类型，对于l1和l2点，线性化的平面内和平面外频率之间的差异有可能以线性化的形式呈现halo/lissayjous轨道，因此对地月共线平动点附近的五种类型的扰动周期轨道进行动力学行为分析；
[0080]
由于稳定/不稳定流形的存在，不可能定位精确的有界轨道。一个小的偏差可能导致路径在几圈后就离开有界轨道。然而，几圈就足以揭示地月共线振动点轨道的特征
[0081]
采用stk/astrogator五组地月halo/lissayjous扰动轨道的初始条件，并用式(16)进行积分，积分时长为路径的3-4圈。历元时刻为2022年10月1日，12：00：00.000图4-8中实线为halo/lissayjous轨道在x-y平面上的投影，红星是轨道与x-y平面的交点，是poincar
é
截面上的点。
[0082]
在lisc中，halo扰动轨道的poincar
é
截面，如图4(b)、6(b)和8(b)所示。仍然与图3(a)所示的典型周期轨道类型相似，虽然一个截面上的点不重叠，但它们仍然彼此紧密相连，并且围绕yi轴对称分布，表现出近似的周期特征。然而，在grc中，相同的halo扰动轨道的poincar
é
截面如图4(a)、6(a)和8(a)所示，失去了周期性的特征。造成grc和lisc中巨大不同现象的主要因素是平动点的运动。在lisc中，坐标系的原点是平动点，它们不再是惯性点，而是随月球半径变化的时间变化点。
[0083]
由此，基于poincar
é
截面分析，新的拼接点由继承的周期部分和扰动部分两部分组成，继承的周期部分可以采用在crtbp模型中的轨道。当轨道在lisc中表示时，只有周期部分显示，因为lisc的起源代表了平动点的运动。扰动部分可以通过求解grc中的li(i＝1,2,3)点的瞬时状态来实现。li(i＝1,2,3))点的瞬时状态包含了由月球偏心度和来自太阳系的引力影响引起的扰动。上述分析将复杂的确定初始条件问题转化为获得crtbp模型中li(i＝1,2,3)的瞬时状态和基本初始条件。因此，由于得到了li(i＝1,2,3)的瞬时状态，可以解决确定初始条件的问题。
[0084]
步骤3:通过poincar
é
截面分析确定平动点瞬时会合坐标系(lisc)下crtbp模型中每段轨迹弧的初始条件
[0085]
由lisc下的crtbp模型生成的轨迹被离散为m个“拼接点”，由m-1个轨迹弧连接。拼接点为一个六维向量xk(k＝1,2,...,m)，每个弧的积分时间为τk(k＝1,2,...,m)，每个点之间的时间间隔为t
k-1
＝τ
k-τ
k-1
，变量向量
[0086]
由于在crtbp模型中，halo/lissayjous轨道在x-y平面上都是周期性的，因此通过在x-y平面上定义的相位角η来定位这些拼接点。η是一个非几何概念下的角度，类似于平近点角的概念。当航天器与x-z平面相交并沿正y方向运动时，相位角为零。相位角沿轨道运动方向增大，其定义为
[0087][0088]
其中p是halo/lissayjous轨道的x-y周期，图9显示了一些不均匀间隔的η值。步骤4:确定平动点的瞬时状态
[0089]
从几何上看，平动点属于月球绕地球的瞬时运动平面，因此到地球和月球的距离满足与crtbp中相同的关系。太阳系天体以及月球绕地球的非圆周运动的影响，使这些点无法成为相对平衡的点。在crtbp中，li(i＝1,2,3)与月球之间的关系参数λi(i＝1,2,3)是仅与质量参数相关的参数。由于月球的偏心度和来自太阳系的引力的影响，λi成为缓慢的时变参数。然而，λi相对于初始时期和时间跨度的变化非常小，这仍然可以被认为是每个短轨迹弧的一个常数。同时，假设li(i＝1,2,3)位置与月球之间的关系为
[0090]rli
＝λrmꢀꢀꢀ
(19)
[0091]
其中rm是在grc中从地球中心到月球中心的位置向量
[0092]
由于平动点是运动方程如等式(17)的稳定解，根据式(19)，其解符合grc中的r
li
＝[λ
irm
(t0),0,0]和采用牛顿迭代过程来确定每个短轨迹弧的λi值,从crtbp中选取初始猜测值λ
i0
。对于第k个弧，由于平动点的初始猜测解不够精确,它们在积分一个时间间隔τk(k＝1,2,...,m)后不能保持平稳，p表示为积分后的解。应使得式(20)的值最小以此确保解符合[λ
irm
(t0),0,0]。
[0093][0094]
其中，t
k-1
为第k个弧的初始时间，τk为第k个弧的时间间隔，可以得到新的λ
i,j+1
。
[0095][0096]
其中δλ
i,j
迭代最小二乘法中的一个小值，超过6次的迭代可以提供理想的结果。步骤5:基于平动点的瞬时状态确定地心会合坐标系(grc)下每段轨迹弧的初始条件,连接拼接点
[0097]
从步骤3中选取历元时刻和拼接点后，根据lisc与grc之间对应的时间依赖关系，将每个拼接点转换为grc中对应的点，具有li(i＝1,2,3)点固定的crtbp模型中，将li(i＝1,2,3)的瞬时状态与继承的初始条件叠加，可以得到更为精确的初始条件。变换关系如图10所示。
[0098]rk
＝r
li
+rkꢀꢀꢀ
(22)
[0099]rk
表示grc中第k个拼接点的位置向量，r
li
表示grc中从地球中心到li(i＝1,2,3)的位置向量，rk表示grc中从li(i＝1,2,3)到探测器的位置向量
[0100][0101]
其中是步骤3中拼接点的信息，ω
cg
是grc相对于lisc的角速度向量,显然，根据这两个坐标系的定义有
[0102]
ω
cg
＝0
ꢀꢀꢀ
(24)
[0103]
将式(24)和(19)代入式(23)中，并考虑下式
[0104]
[0105]
其中是从j2000到grc的变换矩阵，并且可以用式(9)表式，可以得到以下关系
[0106][0107]
其中rm和分别表示月球在j2000中的位置和速度矢量，由于两个向量的点积与坐标系无关，因此rm和的点积可以根据j2000中的星历数据直接计算出来。连接所有拼接点，假设沿每个段的时间不是常数的，约束应该满足新的状态集φk(xk)＝x
k+1
k＝1,2,...,m-1，对于大多数任务来说，最后时刻为τm，因此约束向量为
[0108][0109]
如图10所示，采用序列二次规划法得到了所有内部拼接点的位置和速度都是连续的解。
[0110]
步骤6:采用序列二次规划法得到所有内部拼接点位置速度的解并进行数值模拟
[0111]
首先，在l1/l2点情况的lisc的crtbp模型中生成halo/lissayjous轨道，以及l3sc中l3点的halo轨道，如图12
–
14所示。
[0112]
为了明确地描述拼接点，选择每一圈的相位角η为45
°
，135
°
，225
°
，或315
°
的状态xi作为拼接点的初始条件。初始状态和最终状态也被选择为拼接点。因此，整个轨迹被离散成30个“拼接点”，由29个轨迹弧连接。由于l3点情况下的x-y面频率只有l2点情况下的x-y频率的一半，所以在相同的时间间隔内，只选择了14个“拼接点”。由此确定crtbp模型中所有拼接点的初始条件。
[0113]
其次，基于步骤4中提出的方法，得到了所有五种情况下的轨迹弧对应的每个λi值。初始猜测值λ
i0
根据式(21)从crtbp中选定为λ
10
＝0.849065766935798，λ
20
＝1.16783268238542和λ
30
＝0.99291206846683,根据式(26)，考虑月球的真实运动，得到新拼接点。因此，根据式(19)通过λ
i0
值确定随时间变化li(i＝1,2,3)点位置。
[0114]
最后，利用序列二次规划得到了所有内部拼接点的位置和速度都是连续的解，如图15-19为五种类型的halo/lissayjous轨道。
[0115]
以上所述为本发明的较佳实施例而已，本发明不应该局限于该实施例和附图所公开的内容。凡是不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。