CN113071712B - 一种月地转移入射变轨策略快速计算方法 - Google Patents

一种月地转移入射变轨策略快速计算方法 Download PDF

Info

Publication number
CN113071712B
CN113071712B CN202110261813.6A CN202110261813A CN113071712B CN 113071712 B CN113071712 B CN 113071712B CN 202110261813 A CN202110261813 A CN 202110261813A CN 113071712 B CN113071712 B CN 113071712B
Authority
CN
China
Prior art keywords
vector
orbit
transfer
orbital
monthly
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110261813.6A
Other languages
English (en)
Other versions
CN113071712A (zh
Inventor
孟占峰
高珊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Spacecraft System Engineering
Original Assignee
Beijing Institute of Spacecraft System Engineering
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Spacecraft System Engineering filed Critical Beijing Institute of Spacecraft System Engineering
Priority to CN202110261813.6A priority Critical patent/CN113071712B/zh
Publication of CN113071712A publication Critical patent/CN113071712A/zh
Application granted granted Critical
Publication of CN113071712B publication Critical patent/CN113071712B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/242—Orbits and trajectories

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供一种月地转移入射变轨策略快速计算方法，考虑了月地转移轨道从交会对接后的环月轨道出发限制，采用无需迭代求解的圆轨道解析方法快速求解得到变轨速度增量Δv1的初值，计算效率高；此外，本发明还首次推导了基于有限推力的状态变量I0和再入点终端瞄准状态变量Ef的解析状态之间转移矩阵的传递关系，也即基于偏导数敏度传递矩阵K得到的修正量ΔI0对状态变量I0进行修正，以此获得满足工程约束的变轨策略高精度数值解——变轨点位置矢量的方位角、变轨点位置矢量的高低角以及变轨发动机工作时长，填补了月地转移入射策略的技术空缺。

Description

一种月地转移入射变轨策略快速计算方法
技术领域
本发明属于深空探测轨道设计领域，尤其涉及一种月地转移入射变轨策略快速计算方法。
背景技术
从环月轨道出发返回地球的月地转移轨道是月球采样返回任务所特有的飞行过程，我国前期探月任务均未涉及，世界上只有美国和苏联完成过类似的飞行过程。中国的月球无人采样返回任务采用了月球轨道交会对接后环月等待，待入射能力最优机会出现时进行月地转移入射的策略，这与美国的载人登月任意时间任意地点返回的方案和苏联的起飞后不进入环月轨道直接返回地球都有很大不同。因此，亟需一种月地转移入射变轨策略快速计算方法，以解决我国月球无人采样返回任务月地转移入射轨道的需求。
发明内容
为解决上述问题，本发明提供一种月地转移入射变轨策略快速计算方法，以满足我国月球无人采样返回任务的变轨需求。
一种月地转移入射变轨策略快速计算方法，包括以下步骤：
S1：采用圆轨道解析方法获取从环月轨道非共面入射到月地转移双曲线轨道的变轨速度增量Δv1；
S2：构建基于有限推力的状态变量I0＝(A,E,Δt)和再入点终端瞄准变量
Figure BDA0002970349020000011
其中，A为变轨点位置矢量的方位角，E为变轨点位置矢量的高低角，Δt为变轨发动机工作时长，Hf为再入点的高度，
Figure BDA0002970349020000021
为再入点的轨道倾角，
Figure BDA0002970349020000022
为再入角；
S3：根据变轨速度增量Δv1获取状态变量I0＝(A,E,Δt)的初始值
Figure BDA0002970349020000023
再根据初始值
Figure BDA0002970349020000024
通过轨道预报得到再入点终端瞄准变量的轨道预报值；
S4：获取再入点终端瞄准变量的轨道预报值与期望值之间的偏差ΔEf，判断偏差ΔEf是否小于设定阈值，若小于，则当前对应的状态变量I0为最终的月地转移入射变轨策略，若不小于，则进入步骤S5；
S5：根据偏差ΔEf获取状态变量I0的修正量ΔI0，其中：
ΔI0＝K+ΔEf
其中，K+为偏导数敏度传递矩阵K的广义逆矩阵，且
Figure BDA0002970349020000025
S6：采用修正量ΔI0对初始值
Figure BDA0002970349020000026
进行修正，再根据得到的修正值重新获取偏差ΔEf与偏导数敏度传递矩阵K，并重复步骤S4～S5，直到偏差ΔEf小于设定阈值。
进一步地，所述初始值
Figure BDA0002970349020000027
的获取方法具体为：
S31：获取变轨速度增量Δv1投影在惯性坐标系下三个坐标轴的分量：
Figure BDA0002970349020000028
S32：根据三轴分量获取初始值
Figure BDA0002970349020000029
Figure BDA00029703490200000210
Figure BDA00029703490200000211
Figure BDA0002970349020000031
其中，m0为变轨前的初始质量，
Figure BDA0002970349020000032
为发动机的秒流量，Isp为发动机的比冲，g为重力加速度常数，exp表示e指数函数，Δm为变轨过程消耗的燃料质量，F为发动机的推力大小。
进一步地，所述变轨速度增量Δv1的获取方法为：
S11：获取初始环月轨道的任一时刻的位置速度矢量(r0,v0)、目标月地转移双曲线轨道对应的
Figure BDA0002970349020000033
矢量，其中，
Figure BDA0002970349020000034
矢量为探测器在目标月地转移双曲线轨道上运行至无穷远处时的剩余速度矢量；
S12：计算初始环月轨道的角动量单位矢量
Figure BDA0002970349020000035
和月地转移双曲线轨道的角动量单位矢量
Figure BDA0002970349020000036
Figure BDA0002970349020000037
Figure BDA0002970349020000038
其中，
Figure BDA0002970349020000039
S13：计算
Figure BDA00029703490200000310
矢量与初始轨道面的夹角δrel：
Figure BDA00029703490200000311
S14：根据初始环月轨道的位置速度矢量(r0,v0)计算得到初始环月轨道的半长轴a，进而得到轨道平均速度
Figure BDA00029703490200000312
其中：
Figure BDA00029703490200000313
其中，μmoon为月球引力常数；
S15：计算归一化的速度大小
Figure BDA0002970349020000041
并构造如下辅助变量：
Figure BDA0002970349020000042
同时，建立表征
Figure BDA0002970349020000043
矢量的速度大小K和方向A的函数w(K,A)：
Figure BDA0002970349020000044
S16：计算月地转移双曲线轨道的飞行速度矢量的反方向
Figure BDA0002970349020000045
与变轨点位置矢量r1的夹角β：
Figure BDA0002970349020000046
同时，构造如下辅助变量：
Figure BDA0002970349020000047
Figure BDA0002970349020000048
e2＝(K2+1)2-K2(K2+2)cos2β
其中，ψ表示
Figure BDA0002970349020000049
矢量与变轨点位置矢量r1的夹角，β表示月地转移双曲线轨道的飞行速度矢量
Figure BDA00029703490200000410
与变轨点位置矢量r1的夹角，e为月地转移双曲线轨道的偏心率，α表示变轨点位置矢量r1与轨道面交线矢量间的夹角，且轨道面交线矢量为初始环月轨道面与月地转移双曲线轨道面的交线所在的矢量；
S17：构造轨道面交线矢量方向上的单位矢量
Figure BDA00029703490200000411
Figure BDA0002970349020000051
S18：将
Figure BDA0002970349020000052
矢量绕单位矢量
Figure BDA0002970349020000053
旋转δrel角度，再绕角动量单位矢量
Figure BDA0002970349020000054
旋转-α角度后，得到变轨点位置矢量r1方向上的单位矢量
Figure BDA0002970349020000055
Figure BDA0002970349020000056
S19：采用初始环月轨道的半长轴近似表示变轨点位置矢量r1的大小如下：
|r1|＝a
进而得到变轨点位置矢量
Figure BDA0002970349020000057
S20：将单位矢量
Figure BDA0002970349020000058
绕角动量单位矢量
Figure BDA0002970349020000059
旋转π/2角度，得到初始环月轨道的的飞行速度单位矢量
Figure BDA00029703490200000510
Figure BDA00029703490200000511
S21：根据飞行速度单位矢量
Figure BDA00029703490200000512
的大小得到初始环月轨道的飞行速度矢量
Figure BDA00029703490200000513
Figure BDA00029703490200000514
Figure BDA00029703490200000515
S22：将单位矢量
Figure BDA00029703490200000516
绕角动量单位矢量
Figure BDA00029703490200000517
轴旋转β角后，得到月地转移双曲线轨道的飞行速度单位矢量
Figure BDA00029703490200000518
Figure BDA00029703490200000519
同时，月地转移双曲线轨道的飞行速度矢量
Figure BDA00029703490200000520
的大小为：
Figure BDA00029703490200000521
则得到月地转移双曲线轨道的飞行速度矢量
Figure BDA0002970349020000061
Figure BDA0002970349020000062
S23：根据飞行速度矢量
Figure BDA0002970349020000063
与飞行速度矢量
Figure BDA0002970349020000064
得到变轨速度增量Δv1：
Figure BDA0002970349020000065
进一步地，根据链式求导传递法则，得到所述偏导数敏度传递矩阵K的表达式：
Figure BDA0002970349020000066
其中，
Figure BDA0002970349020000067
为月地转移双曲线轨道终端地球固连系下的位置速度，
Figure BDA0002970349020000068
为月地转移双曲线轨道惯性系下的位置速度，
Figure BDA0002970349020000069
为发动机关机点地心惯性系下的位置速度，
Figure BDA00029703490200000610
为发动机关机点月心惯性系下的位置速度。
进一步地，步骤S4中所述的设定阈值为10－3。
有益效果：
本发明提供一种月地转移入射变轨策略快速计算方法，考虑了月地转移轨道从交会对接后的环月轨道出发限制，采用无需迭代求解的圆轨道解析方法快速求解得到变轨速度增量Δv1的初值，计算效率高；此外，本发明还首次推导了基于有限推力的状态变量I0和再入点终端瞄准状态变量Ef的解析状态之间转移矩阵的传递关系，也即基于偏导数敏度传递矩阵K得到的修正量ΔI0对状态变量I0进行修正，以此获得满足工程约束的变轨策略高精度数值解——变轨点位置矢量的方位角、变轨点位置矢量的高低角以及变轨发动机工作时长，填补了月地转移入射策略的技术空缺。
附图说明
图1为本发明提供的变轨示意图。
具体实施方式
为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述。
本发明以我国月球无人采样返回任务月地转移入射轨道策略设计为需求，提出了一种月地转移入射变轨策略快速计算方法，该方法采用解析算法可以直接求解得到从环月轨道非共面入射进入月地转移轨道的能量最优变轨点位置和变轨速度增量大小；然后采用全解析状态转移矩阵，建立了有限推力高精度模型下，设计变量与目标变量的偏差解析传递关系；以解析初值为基础，采用高精度数值积分模型利用所建立的偏差传递关系就可以得到真实引力场下的高精度能量最优变轨结果。
具体的，一种月地转移入射变轨策略快速计算方法，其特征在于，包括以下步骤：
S1：采用圆轨道解析方法获取从环月轨道非共面入射到月地转移双曲线轨道的变轨速度增量Δv1。
S2：构建状态变量I0＝(A,E,Δt)和再入点终端瞄准变量
Figure BDA0002970349020000071
其中，A为变轨点位置矢量在J2000惯性标系下的方位角，E为变轨点位置矢量在J2000惯性标系下的高低角，Δt为变轨发动机工作时长，Hf为再入点的高度，
Figure BDA0002970349020000072
为再入点的轨道倾角，
Figure BDA0002970349020000073
为再入角。
S3：根据变轨速度增量Δv1获取状态变量I0＝(A,E,Δt)的初始值
Figure BDA0002970349020000074
再根据初始值
Figure BDA0002970349020000075
通过轨道预报得到再入点终端瞄准变量的轨道预报值；
具体的，初始值
Figure BDA0002970349020000081
的获取方法具体为：
S31：获取变轨速度增量Δv1投影在惯性坐标系下三个坐标轴的分量：
Figure BDA0002970349020000082
其中，T表示转置；
S32：根据三轴分量获取初始值
Figure BDA0002970349020000083
Figure BDA0002970349020000084
Figure BDA0002970349020000085
Figure BDA0002970349020000086
其中，上标IG是Initial Guess(初值猜测)的缩写，m0为变轨前的初始质量，
Figure BDA0002970349020000087
为发动机的秒流量，Isp为发动机的比冲，g为重力加速度常数，g＝9.80665m/s2，exp表示e指数函数，Δm为变轨过程消耗的燃料质量，F为发动机的推力大小。
S4：获取再入点终端瞄准变量的轨道预报值与期望值之间的偏差ΔEf，判断偏差ΔEf是否小于设定阈值，如10－3，若小于，则当前对应的状态变量I0为最终的月地转移入射变轨策略，若不小于，则进入步骤S5。
S5：根据偏差ΔEf获取状态变量I0的修正量ΔI0，其中：
ΔI0＝K+ΔEf
其中，K+为偏导数敏度传递矩阵K的广义逆矩阵，且
Figure BDA0002970349020000088
S6：采用修正量ΔI0对初始值
Figure BDA0002970349020000091
进行修正，如两者相加，再根据得到的修正值重新获取偏差ΔEf与偏导数敏度传递矩阵K，并重复步骤S4～S5，直到偏差ΔEf小于设定阈值。
进一步地，变轨速度增量Δv1的获取方法为：
S11：获取初始环月轨道的任一时刻的位置速度矢量(r0,v0)、目标月地转移双曲线轨道对应的
Figure BDA0002970349020000092
矢量，其中，如图1所示，
Figure BDA0002970349020000093
矢量为探测器在目标月地转移双曲线轨道上运行至无穷远处时的剩余速度矢量；
S12：计算初始环月轨道的角动量单位矢量
Figure BDA0002970349020000094
和月地转移双曲线轨道的角动量单位矢量
Figure BDA0002970349020000095
Figure BDA0002970349020000096
Figure BDA0002970349020000097
其中，
Figure BDA0002970349020000098
S13：计算
Figure BDA0002970349020000099
矢量与初始轨道面的夹角δrel：
Figure BDA00029703490200000910
S14：根据初始环月轨道的位置速度矢量(r0,v0)计算得到初始环月轨道的半长轴a，进而得到轨道平均速度
Figure BDA00029703490200000913
其中：
Figure BDA00029703490200000911
其中，μmoon为月球引力常数；
S15：计算归一化的速度大小
Figure BDA00029703490200000912
并构造如下辅助变量：
Figure BDA0002970349020000101
同时，建立表征
Figure BDA0002970349020000102
矢量的速度大小K和方向A的函数w(K,A)：
Figure BDA0002970349020000103
S16：计算月地转移双曲线轨道的飞行速度矢量的反方向
Figure BDA0002970349020000104
与变轨点位置矢量r1的夹角β：
Figure BDA0002970349020000105
同时，构造如下辅助变量：
Figure BDA0002970349020000106
Figure BDA0002970349020000107
e2＝(K2+1)2-K2(K2+2)cos2β
其中，ψ表示
Figure BDA0002970349020000108
矢量与变轨点位置矢量r1的夹角，β表示月地转移双曲线轨道的飞行速度矢量
Figure BDA0002970349020000109
与变轨点位置矢量r1的夹角，e为月地转移双曲线轨道的偏心率，α表示变轨点位置矢量r1与轨道面交线矢量间的夹角，且轨道面交线矢量为初始环月轨道面与月地转移双曲线轨道面的交线所在的矢量；
下面根据前述计算得到的标量和角度值，转换为位置和速度的矢量表达形式，计算变轨点位置矢量r1、变轨前速度矢量
Figure BDA00029703490200001010
以及变轨后速度矢量
Figure BDA00029703490200001011
具体参见如下步骤S17～S23。
S17：构造轨道面交线矢量方向上的单位矢量
Figure BDA0002970349020000111
Figure BDA0002970349020000112
S18：将
Figure BDA0002970349020000113
矢量绕单位矢量
Figure BDA0002970349020000114
旋转δrel角度，再绕角动量单位矢量
Figure BDA0002970349020000115
旋转-α角度后，得到变轨点位置矢量r1方向上的单位矢量
Figure BDA0002970349020000116
Figure BDA0002970349020000117
S19：采用初始环月轨道的半长轴近似表示变轨点位置矢量r1的大小如下：
|r1|＝a
进而得到变轨点位置矢量
Figure BDA0002970349020000118
S20：将单位矢量
Figure BDA0002970349020000119
绕角动量单位矢量
Figure BDA00029703490200001110
旋转π/2角度，得到初始环月轨道的的飞行速度单位矢量
Figure BDA00029703490200001111
Figure BDA00029703490200001112
S21：根据飞行速度单位矢量
Figure BDA00029703490200001113
的大小得到初始环月轨道的飞行速度矢量
Figure BDA00029703490200001114
Figure BDA00029703490200001115
Figure BDA00029703490200001116
S22：将单位矢量
Figure BDA00029703490200001117
绕角动量单位矢量
Figure BDA00029703490200001118
轴旋转β角后，得到月地转移双曲线轨道的飞行速度单位矢量
Figure BDA00029703490200001119
Figure BDA00029703490200001120
同时，月地转移双曲线轨道的飞行速度矢量
Figure BDA0002970349020000121
的大小为：
Figure BDA0002970349020000122
则得到月地转移双曲线轨道的飞行速度矢量
Figure BDA0002970349020000123
Figure BDA0002970349020000124
S23：根据飞行速度矢量
Figure BDA0002970349020000125
与飞行速度矢量
Figure BDA0002970349020000126
得到变轨速度增量Δv1：
Figure BDA0002970349020000127
需要说明的是，在得到变轨速度增量Δv1后，本发明再以变轨速度增量Δv1为初值，计算实际椭圆轨道在有限推力和复杂摄动引力场下的最优变轨点位置、发动机的推力方向、发动机的推力时长等月地转移入射变轨策略。
进一步地，本发明将脉冲变轨模型转化为有限推力变轨模型，需要建立有限推力设计状态变量I0＝(A,E,Δt)和再入点终端瞄准变量
Figure BDA0002970349020000128
的偏导数敏度传递矩阵
Figure BDA0002970349020000129
表达式。
初始状态变量偏差ΔI0和再入点终端瞄准变量偏差ΔEf满足ΔEf＝KΔI0，则可以推导出ΔI0＝K+ΔEf；同时，根据链式求导传递法则，偏导数敏度传递矩阵
Figure BDA00029703490200001210
表达式如下：
Figure BDA00029703490200001211
其中，
Figure BDA00029703490200001212
为月地转移双曲线轨道终端地球固连系下的位置速度，
Figure BDA00029703490200001213
为月地转移双曲线轨道惯性系下的位置速度，
Figure BDA00029703490200001214
为发动机关机点地心惯性系下的位置速度，
Figure BDA00029703490200001215
为发动机关机点月心惯性系下的位置速度。
下面对偏导数敏度传递矩阵
Figure BDA00029703490200001216
表达式中的各项计算步骤分别进行说明。
(2.1)计算变轨速度增量Δv1对有限推力参数I0＝(A,E,Δt)的偏导数关系
Figure BDA0002970349020000131
Figure BDA0002970349020000132
其中，
Figure BDA0002970349020000133
(2.2)计算发动机关机点位置速度
Figure BDA0002970349020000134
对变轨速度增量Δv1的偏导数关系
Figure BDA0002970349020000135
Figure BDA0002970349020000136
其中，m(t)表示变轨过程中探测器质量关于时间的变化函数；
M(t)＝m(t)/m0为质量剩余比例，为归一化函数，表示变轨过程中质量相比初始质量的比值，0≤M(t)≤1；
Mf表示变轨结束时刻的质量剩余比例，即：Mf＝M(tf)；
B＝m0E/lnMf；
E表示Δv1投影坐标系相对惯性坐标系的转换矩阵，本发明Δv1在惯性坐标系下投影，因此E＝I3×3，为3×3单位矢量矩阵；
(2.3)计算终端惯性系位置速度
Figure BDA0002970349020000137
相对初始位置速度
Figure BDA0002970349020000138
的偏导数关系
Figure BDA0002970349020000139
该数值可以由
Figure BDA00029703490200001310
到
Figure BDA00029703490200001311
轨道预报计算获得。
(2.4)计算月地转移双曲线轨道终端固连系状态
Figure BDA00029703490200001312
相对终端惯性系位置速度
Figure BDA00029703490200001313
的传递关系
Figure BDA00029703490200001314
Figure BDA0002970349020000141
其中，
Figure BDA0002970349020000142
为地心惯性坐标系到地心固连坐标系的坐标转换矩阵；
Figure BDA0002970349020000143
为地心惯性坐标系到地心固连坐标系的坐标转换矩阵变化率矩阵。
(2.5)计算月地转移双曲线轨道终端状态，即月地转移双曲线轨道的终点处的
Figure BDA0002970349020000144
对终端位置速度
Figure BDA0002970349020000145
的偏导数关系
Figure BDA0002970349020000146
Figure BDA0002970349020000147
其中，
Figure BDA0002970349020000148
Figure BDA0002970349020000149
Figure BDA00029703490200001410
各项的具体表达为：
Figure BDA00029703490200001411
Figure BDA0002970349020000151
将(2.1)～(2.5)的计算结果带入偏导数敏度传递矩阵
Figure BDA0002970349020000152
的表达式，就可以计算出微分修正敏度矩阵
Figure BDA0002970349020000153
的结果。
综上，本发明采取的技术解决方案如下：
首先获取探测器的环月初始轨道位置速度矢量(r0,v0)、目标月地转移轨道的
Figure BDA0002970349020000154
矢量和再入点目标参数(高度、倾角和再入角)；其次，根据初始环月轨道参数和目标月地转移轨道的
Figure BDA0002970349020000155
矢量，利用圆轨道解析方法，基于变轨点位置计算变轨速度增量初值；然后，以变轨速度增量初值，计算实际椭圆轨道在有限推力和复杂摄动引力场下的最优变轨点位置和变轨过程。
与现有技术相比，本发明的优点是：1)方法考虑了月地转移轨道从交会对接后的环月轨道出发限制；2)采用无需迭代求解的解析方法快速求解得到初值，计算效率高；3)首次推导了有限推力和再入点终端瞄准状态的解析状态转移矩阵的传递关系，可以获得满足工程约束的高精度数值解。
当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当然可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (4)

1.一种月地转移入射变轨策略快速计算方法，其特征在于，包括以下步骤：
S1：采用圆轨道解析方法获取从环月轨道非共面入射到月地转移双曲线轨道的变轨速度增量Δv1；
S2：构建基于有限推力的状态变量I0＝(A,E,Δt)和再入点终端瞄准变量
Figure FDA0003587412330000011
其中，A为变轨点位置矢量的方位角，E为变轨点位置矢量的高低角，Δt为变轨发动机工作时长，Hf为再入点的高度，
Figure FDA0003587412330000012
为再入点的轨道倾角，
Figure FDA0003587412330000013
为再入角；
S3：根据变轨速度增量Δv1获取状态变量I0＝(A,E,Δt)的初始值
Figure FDA0003587412330000014
再根据初始值
Figure FDA0003587412330000015
通过轨道预报得到再入点终端瞄准变量的轨道预报值；所述初始值
Figure FDA0003587412330000016
的获取方法具体为：
S31：获取变轨速度增量Δv1投影在惯性坐标系下三个坐标轴的分量：
Figure FDA0003587412330000017
S32：根据三轴分量获取初始值
Figure FDA0003587412330000018
Figure FDA0003587412330000019
Figure FDA00035874123300000110
Figure FDA00035874123300000111
其中，m0为变轨前的初始质量，
Figure FDA00035874123300000112
为发动机的秒流量，Isp为发动机的比冲，g为重力加速度常数，exp表示e指数函数，Δm为变轨过程消耗的燃料质量，F为发动机的推力大小；
S4：获取再入点终端瞄准变量的轨道预报值与期望值之间的偏差ΔEf，判断偏差ΔEf是否小于设定阈值，若小于，则当前对应的状态变量I0为最终的月地转移入射变轨策略，若不小于，则进入步骤S5；
S5：根据偏差ΔEf获取状态变量I0的修正量ΔI0，其中：
ΔI0＝K+ΔEf
其中，K+为偏导数敏度传递矩阵K的广义逆矩阵，且
Figure FDA0003587412330000021
S6：采用修正量ΔI0对初始值
Figure FDA0003587412330000022
进行修正，再根据得到的修正值重新获取偏差ΔEf与偏导数敏度传递矩阵K，并重复步骤S4～S5，直到偏差ΔEf小于设定阈值。
2.如权利要求1所述的一种月地转移入射变轨策略快速计算方法，其特征在于，所述变轨速度增量Δv1的获取方法为：
S11：获取初始环月轨道的任一时刻的位置速度矢量(r0,v0)、目标月地转移双曲线轨道对应的
Figure FDA0003587412330000023
矢量，其中，
Figure FDA0003587412330000024
矢量为探测器在目标月地转移双曲线轨道上运行至无穷远处时的剩余速度矢量；
S12：计算初始环月轨道的角动量单位矢量
Figure FDA0003587412330000025
和月地转移双曲线轨道的角动量单位矢量
Figure FDA0003587412330000026
Figure FDA0003587412330000027
Figure FDA0003587412330000028
其中，
Figure FDA0003587412330000029
S13：计算
Figure FDA00035874123300000210
矢量与初始轨道面的夹角δrel：
Figure FDA0003587412330000031
S14：根据初始环月轨道的位置速度矢量(r0,v0)计算得到初始环月轨道的半长轴a，进而得到轨道平均速度
Figure FDA0003587412330000032
其中：
Figure FDA0003587412330000033
其中，μmoon为月球引力常数；
S15：计算归一化的速度大小
Figure FDA0003587412330000034
并构造如下辅助变量：
Figure FDA0003587412330000035
同时，建立表征
Figure FDA0003587412330000036
矢量的速度大小K和方向A的函数w(K,A)：
Figure FDA0003587412330000037
S16：计算月地转移双曲线轨道的飞行速度矢量的反方向
Figure FDA0003587412330000038
与变轨点位置矢量r1的夹角β：
Figure FDA0003587412330000039
同时，构造如下辅助变量：
Figure FDA00035874123300000310
Figure FDA00035874123300000311
e2＝(K2+1)2-K2(K2+2)cos2β
其中，ψ表示
Figure FDA0003587412330000041
矢量与变轨点位置矢量r1的夹角，β表示月地转移双曲线轨道的飞行速度矢量
Figure FDA0003587412330000042
与变轨点位置矢量r1的夹角，e为月地转移双曲线轨道的偏心率，α表示变轨点位置矢量r1与轨道面交线矢量间的夹角，且轨道面交线矢量为初始环月轨道面与月地转移双曲线轨道面的交线所在的矢量；
S17：构造轨道面交线矢量方向上的单位矢量
Figure FDA0003587412330000043
Figure FDA0003587412330000044
S18：将
Figure FDA0003587412330000045
矢量绕单位矢量
Figure FDA0003587412330000046
旋转δrel角度，再绕角动量单位矢量
Figure FDA0003587412330000047
旋转-α角度后，得到变轨点位置矢量r1方向上的单位矢量
Figure FDA0003587412330000048
Figure FDA0003587412330000049
S19：采用初始环月轨道的半长轴近似表示变轨点位置矢量r1的大小如下：
|r1|＝a
进而得到变轨点位置矢量
Figure FDA00035874123300000410
S20：将单位矢量
Figure FDA00035874123300000411
绕角动量单位矢量
Figure FDA00035874123300000412
旋转π/2角度，得到初始环月轨道的飞行速度单位矢量
Figure FDA00035874123300000413
Figure FDA00035874123300000414
S21：根据飞行速度单位矢量
Figure FDA00035874123300000415
的大小得到初始环月轨道的飞行速度矢量
Figure FDA00035874123300000416
Figure FDA00035874123300000417
Figure FDA00035874123300000418
S22：将单位矢量
Figure FDA0003587412330000051
绕角动量单位矢量
Figure FDA0003587412330000052
轴旋转β角后，得到月地转移双曲线轨道的飞行速度单位矢量
Figure FDA0003587412330000053
Figure FDA0003587412330000054
同时，月地转移双曲线轨道的飞行速度矢量
Figure FDA0003587412330000055
的大小为：
Figure FDA0003587412330000056
则得到月地转移双曲线轨道的飞行速度矢量
Figure FDA0003587412330000057
Figure FDA0003587412330000058
S23：根据飞行速度矢量
Figure FDA0003587412330000059
与飞行速度矢量
Figure FDA00035874123300000510
得到变轨速度增量Δv1：
Figure FDA00035874123300000511
3.如权利要求1所述的一种月地转移入射变轨策略快速计算方法，其特征在于，根据链式求导传递法则，得到所述偏导数敏度传递矩阵K的表达式：
Figure FDA00035874123300000512
其中，
Figure FDA00035874123300000513
为月地转移双曲线轨道终端地球固连系下的位置速度，
Figure FDA00035874123300000514
为月地转移双曲线轨道惯性系下的位置速度，
Figure FDA00035874123300000515
为发动机关机点地心惯性系下的位置速度，
Figure FDA00035874123300000516
为发动机关机点月心惯性系下的位置速度。
4.如权利要求1所述的一种月地转移入射变轨策略快速计算方法，其特征在于，步骤S4中所述的设定阈值为10－3。
CN202110261813.6A 2021-03-10 2021-03-10 一种月地转移入射变轨策略快速计算方法 Active CN113071712B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110261813.6A CN113071712B (zh) 2021-03-10 2021-03-10 一种月地转移入射变轨策略快速计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110261813.6A CN113071712B (zh) 2021-03-10 2021-03-10 一种月地转移入射变轨策略快速计算方法

Publications (2)

Publication Number Publication Date
CN113071712A CN113071712A (zh) 2021-07-06
CN113071712B true CN113071712B (zh) 2022-07-29

Family

ID=76612729

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110261813.6A Active CN113071712B (zh) 2021-03-10 2021-03-10 一种月地转移入射变轨策略快速计算方法

Country Status (1)

Country Link
CN (1) CN113071712B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109625323A (zh) * 2018-11-09 2019-04-16 中国科学院空间应用工程与技术中心 一种卫星化学推进变轨方法及系统
CN110015445A (zh) * 2019-02-15 2019-07-16 北京空间飞行器总体设计部 一种地月L2点Halo轨道维持方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
AU7357198A (en) * 1997-03-25 1998-10-20 Edward A. Belbruno Low energy method for changing the inclinations of orbiting satellites using weak stability boundaries and a computer process for implementing same
US6237876B1 (en) * 2000-07-28 2001-05-29 Space Systems/Loral, Inc. Methods for using satellite state vector prediction to provide three-axis satellite attitude control
CN101354251B (zh) * 2008-09-12 2010-09-15 航天东方红卫星有限公司 一种深空探测器等效转移轨道确定方法
CN110096726B (zh) * 2019-02-21 2023-08-11 上海卫星工程研究所 基于月球借力的geo卫星应急转移轨道快速优化设计方法
CN110733667B (zh) * 2019-09-29 2021-06-18 北京空间技术研制试验中心 地月平动点轨道间转移设计方法
CN110704952B (zh) * 2019-09-30 2022-09-09 中国人民解放军国防科技大学 一种月地三脉冲返回轨道速度增量分析方法
CN110909461B (zh) * 2019-11-13 2021-09-17 清华大学 基于可达集概念的地月/月地直接转移轨道设计方法
CN111361762B (zh) * 2020-03-04 2020-12-11 北京空间飞行器总体设计部 一种地月转移轨道发动机试喷方法
CN111460614B (zh) * 2020-03-04 2020-12-11 北京空间飞行器总体设计部 一种地月l2点转移轨道中途修正方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109625323A (zh) * 2018-11-09 2019-04-16 中国科学院空间应用工程与技术中心 一种卫星化学推进变轨方法及系统
CN110015445A (zh) * 2019-02-15 2019-07-16 北京空间飞行器总体设计部 一种地月L2点Halo轨道维持方法

Also Published As

Publication number Publication date
CN113071712A (zh) 2021-07-06

Similar Documents

Publication Publication Date Title
CN109375648B (zh) 一种多约束条件下椭圆轨道卫星编队构形初始化方法
Kluever Simple guidance scheme for low-thrust orbit transfers
Huang et al. Post-capture attitude control for a tethered space robot–target combination system
CN106679674B (zh) 基于星历模型的地月L2点Halo轨道阴影分析方法
CN112607065B (zh) 一种基于电推进系统的高精度相位控制方法
Dang et al. Modeling and analysis of relative hovering control for spacecraft
CN111552003B (zh) 基于球卫星编队的小行星引力场全自主测量系统及方法
CN112989496B (zh) 航天器制导方法、装置、电子设备及存储介质
CN109240340B (zh) 一种基于拟周期轨道的洛伦兹力多星编队构型方法
CN113602535A (zh) 一种微纳卫星在轨自主交会控制的方法及计算机设备
CN111338367A (zh) 一种偏心率冻结同轨双脉冲控制的中间轨道确定方法
CN111460614B (zh) 一种地月l2点转移轨道中途修正方法
Sah et al. Design development of debris chaser small satellite with robotic manipulators for debris removal
CN113071712B (zh) 一种月地转移入射变轨策略快速计算方法
Ito et al. Sounding rocket SS-520 as a CubeSat launch vehicle
CN113093246A (zh) 地面多目标点成像快速判定及任务参数计算方法
Lu et al. Decentralized closed-loop optimization for 6-DOF self-assembly satellites
Han et al. Engineering Optimization Method of Orbit Transfer Strategy for All-electric Propulsion Satellites
CN112009727B (zh) 平动点轨道最优小推力转移分段设计方法
Broschart et al. Shadow navigation support at jpl for the rosetta landing on comet 67p/churyumov-gerasimenko
Lu et al. Design and optimization of low-energy transfer orbit to Mars with multi-body environment
CN111475767A (zh) 考虑地球自转影响的最小能量弹道严格构造方法
Lee et al. Manifold-Based Space Mission Design with Poincaré Filtering Algorithm
CN113282097B (zh) 一种geo博弈航天器相对位置非球形摄动误差的控制方法
Zhang et al. Space target surveillance based on non-linear model predictive control

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant