CN109774974B - 一种用于空间碎片接近的轨道设计方法 - Google Patents

一种用于空间碎片接近的轨道设计方法 Download PDF

Info

Publication number
CN109774974B
CN109774974B CN201910096002.8A CN201910096002A CN109774974B CN 109774974 B CN109774974 B CN 109774974B CN 201910096002 A CN201910096002 A CN 201910096002A CN 109774974 B CN109774974 B CN 109774974B
Authority
CN
China
Prior art keywords
intersection
track
point
intersection point
points
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
CN201910096002.8A
Other languages
English (en)
Other versions
CN109774974A (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.)
Shanghai Zhongkechen New Satellite Technology Co ltd
Original Assignee
Shanghai Engineering Center for Microsatellites
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 Shanghai Engineering Center for Microsatellites filed Critical Shanghai Engineering Center for Microsatellites
Priority to CN201910096002.8A priority Critical patent/CN109774974B/zh
Publication of CN109774974A publication Critical patent/CN109774974A/zh
Application granted granted Critical
Publication of CN109774974B publication Critical patent/CN109774974B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明提供了一种用于空间碎片接近的轨道设计方法，包括设定机动卫星的停泊轨道；将所述停泊轨道和所述目标轨道的交线地心距差足够小的交线点取出作为可能交会点，根据所述可能交会点取出对应的可能轨控点；在无时间的约束下，求解每个所述可能交会点和与其相应的所述可能轨控点之间的最小速度增量，从而每个所述可能交会点均可得到与其对应的速度增量最小值和最优轨控点；结合相位调整对速度增量的影响，优选出最佳交会点和最佳轨控点；最后将所述机动卫星与目标碎片在接近时的距离矢量要求叠加到所述最佳交会点上，迭代计算最终所需速度增量，得到最优转移轨道。

Description

一种用于空间碎片接近的轨道设计方法
技术领域
本发明涉及航天技术领域，尤其涉及一种用于空间碎片接近的轨道设计方法。
背景技术
自1957年苏联第一颗卫星升空至2018年1月为止，全球已有5400次的发射，短短六十年，空间技术得到了快速发展，与此同时，太空中也遗留了越来越多的空间碎片，轨道变得越来越拥挤，碰撞事件频发：1991年12月，当时俄罗斯一颗失效卫星“宇宙1934”撞上了本国另一颗卫星“宇宙926”；1996年7月，法国Cerise电子侦察卫星与“Ariane V16”末级火箭残骸相撞；2005年1月，中国“长征4号”火箭碎片和美国发射的“雷神”火箭末级残骸发生碰撞；2009年2月，美国的通信卫星“铱星33”与俄罗斯已经报废的军用卫星“宇宙2251”号碰撞。这些碰撞又产生了数以千计的碎片，给原本恶劣的轨道环境雪上加霜。依据美国空间监测网站(Space Surveillance Network，SSN)的数据，至2018年1月，小于1厘米的碎片超过1.66亿个，1-10厘米直径大小的碎片约有75万个，大于10厘米的碎片约有29000个。
上述碎片的绝大多数都分布在近地轨道上，虽然有一部分碎片可以在大气阻力的作用下自然陨落，但大部分则继续呆在太空中。近地轨道是人类使用最为频繁的轨道，诸如SpaceX、OneWeb等公司都提出了在近地轨道上建立卫星颗数以百计甚至千计的星座计划，这也意味着近地轨道将进一步拥挤。为了能够保证今后近地轨道的合理科学的继续使用，一是令后续发射的卫星在服务结束后能主动离轨，二是采用主动手段移除空间碎片。对于后者，目前已经有用抓钩抓取碎片、使用渔网给碎片增阻、使用大太阳帆等手段，有一些还停留在概念阶段，有一些已经进入太空试验阶段，2018年10月，英国用渔网抓捕小卫星试验取得成功。无论是抓钩、渔网还是太阳帆，都需要主动接近碎片才有效。对于未来专门进行碎片清除任务的卫星来说，很明显，为了提高效率和降低成本，需尽量多的清除碎片，在这也意味如何在有限的燃料下接近更多的碎片。
发明内容
有鉴于上述技术问题，本发明提供了一种用于空间碎片接近的椭圆轨道设计方法，包括：
设定机动卫星的停泊轨道；
将所述停泊轨道和所述目标轨道的交线地心距差足够小的交线点取出作为可能交会点，根据所述可能交会点取出对应的可能轨控点；
在无时间的约束下，求解每个所述可能交会点和与其相应的所述可能轨控点之间的最小速度增量，从而每个所述可能交会点均可得到与其对应的速度增量最小值和最优轨控点；
结合相位调整对速度增量的影响，优选出最佳交会点和最佳轨控点。
将所述机动卫星与目标碎片在接近时的距离矢量要求叠加到所述最佳交会点上，迭代计算最终所需速度增量，得到最优转移轨道。
进一步地，所述停泊轨道被设置为椭圆轨道。
进一步地，所述停泊轨道的参数设置包括：停泊轨道倾角0°、近地点高度400km、远地点高度1400km。
进一步地，所述交线地心距差是停泊轨道交点径向地心距与目标轨道交点径向地心距的差异。
进一步地，将所述交线地心距差小于20km的交线点取出作为所述可能交会点。
进一步地，取出距所述可能交会点最近的一轨数据为所述可能轨控点。
进一步地，从所述可能交会点对应时刻回推一轨，通过插值得到360个所述可能轨控点，每两个所述可能轨控点间隔1度。
进一步地，在求解每个所述可能交会点和与其相应的所述可能轨控点之间的最小速度增量时，通过平面投影法克服所述停泊轨道和所述目标轨道的法向偏差。
进一步地，在求解每个所述可能交会点和与其相应的所述可能轨控点之间的最小速度增量时，所述可能交会点采用相应的虚拟交会点，所述虚拟交会点是事先计算出二体模型与精确模型外推后的偏差矢量，然后将所述偏差矢量叠加到当前交会点后形成的交会点。
本发明的用于空间碎片接近的轨道设计方法，依据近地轨道的摄动特性，设计出相应的椭圆轨道，首先利用摄动，实现停泊轨道与目标轨道两者在轨道面交线方向上距离的接近，然后采用卫星控制，实现机动卫星与目标碎片的接近。通过摄动主控，卫星辅控的方式，只需很少的燃料就能接近碎片，在相同的燃料下能完成更多的接近动作。本发明适用于接近碎片，以及其它对时间要求不高的航天接近任务。
以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。
附图说明
图1是交线运动分解示意图，AB和CD为两个方向的交线地心距差；
图2是本发明的一个实施例中
Figure BDA0001964567660000031
随停泊轨道倾角iP的变化曲线图；
图3是本发明的一个实施例中iP＝0°，iP＝180°两种情况下
Figure BDA0001964567660000032
随目标轨道倾角iT的变化曲线图；
图4是本发明的一个实施例的仿真实验中CZ-2C在交线点时与停泊轨道的交线地心距差的变化曲线图；
图5是本发明的一个实施例的仿真实验中机动卫星与CZ-2C碎片交会接近三维图；
图6是本发明的一个实施例的仿真实验中机动卫星与CZ-2C碎片距离的变化曲线图。
具体实施方式
在本发明的实施方式的描述中，需要理解的是，术语“上”、“下”、“前”、“后”、“左”、“右”、“垂直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对发明的限制。附图为原理图或者概念图，各部分厚度与宽度之间的关系，以及各部分之间的比例关系等等，与其实际值并非完全一致。
本实施例的用于空间碎片接近的椭圆轨道设计方法，将机动卫星的停泊轨道设计成合适的椭圆轨道，在摄动作用下，停泊轨道能与目标轨道的交线地心距差接近甚至相交。将交线地心距差足够小的交线点作为可能交会点，取出距可能交会点最近的一轨数据为可能轨控点。因为交线地心距差足够小，初期计算的时候先设定交线地心距差为0。无时间的约束下，求解可能交会点和可能轨控点的最小速度增量，最终得到与每一个可能交会点对应速度增量最小的最优轨控点。再考虑相位调整对速度增量的影响，最终优选出最佳交会点和最佳轨控点。当机动卫星与目标碎片最终接近距离矢量要求为d，则可将此矢量叠加到最佳交会点上，迭代计算最终所需速度增量，得到最优转移轨道。按照上述的计算方法，最终只需很小的速度增量即可实现碎片接近，具体设计步骤如下：
(一)椭圆停泊轨道设计
1.交线地心距差
机动卫星轨道面与目标轨道面相交的线为交线矢量，交线矢量与停泊轨道面和目标轨道面都有交点，将从节线到交线矢量的角称为交线幅角，目标轨道交线幅角uTc和停泊轨道交线幅角uPc求解如下：
Figure BDA0001964567660000041
Figure BDA0001964567660000042
其中：ΩP为停泊轨道升交点赤经；
iP为停泊轨道倾角；
ΩT为目标轨道升交点赤经；
iT为目标轨道倾角。
停泊轨道交点径向地心距与目标轨道交点径向地心距的差异为交线地心距差，求解如下：
Figure BDA0001964567660000043
其中：
fTc是目标轨道所在交线对应的真近点角值，fTc＝uTc-wT，uTc为目标碎片交线幅角，wT为目标碎片近地点角；
fPc是停泊轨道所在交线对应的真近点角值，fPc＝uPc-wP，uPc为机动卫星交线幅角，wP为机动卫星近地点角；
ΔR0为目标轨道交线幅角uTc，停泊轨道交线幅角uPc时的交线地心距差；
ΔR180为目标轨道交线幅角uTc+180°，停泊轨道交线幅角uPc+180°时的交线地心距差；
aT为目标轨道的半长轴；
aP为停泊轨道的半长轴；
eT为目标轨道的偏心率；
eP为停泊轨道的偏心率。
2.摄动的影响
一阶长期项下，倾角、偏心率、半长轴可认为不变，但升交点赤经会有变化，升交点赤经变化率
Figure BDA0001964567660000044
如下：
Figure BDA0001964567660000045
其中，J2＝1082.63×10-6
Re为地球半径，为6378.14km；
e为对应轨道的偏心率；
a为对应轨道的半长轴；
i为对应轨道的倾角；
n为对应轨道自转角速率，
Figure BDA0001964567660000051
μ＝398600.4418km3/sec2。
近地点角变化率
Figure BDA0001964567660000052
如下：
Figure BDA0001964567660000053
由此可见，由于摄动的作用，机动卫星与目标轨道的升交点赤经会随时间变化：
Figure BDA0001964567660000054
Figure BDA0001964567660000055
其中：ΩP0为停泊轨道初始升交点赤经；
ΩT0为目标轨道初始升交点赤经；
Figure BDA0001964567660000056
为停泊轨道升交点赤经变化率，将停泊轨道参数代入式(4)可得；
Figure BDA0001964567660000057
为目标轨道升交点赤经变化率，将目标轨道参数代入式(4)可得；
Δt为当前时刻距离初始时刻的时间。
由于摄动的关系，停泊轨道和目标轨道近地点角随时间的变化如下：
Figure BDA0001964567660000058
Figure BDA0001964567660000059
其中，wP0为停泊轨道初始近地点角；
Figure BDA00019645676600000510
将停泊轨道参数代入式(5)得到的近地点角变化率；
wT0为目标轨道初始近地点角；
Figure BDA00019645676600000511
将目标轨道参数代入式(5)得到的近地点角变化率。
3.交线真近角变化率和覆盖周期
停泊轨道交线真近点变化率
Figure BDA00019645676600000512
Figure BDA0001964567660000061
其中：
Figure BDA0001964567660000062
Figure BDA0001964567660000063
覆盖周期T为：
Figure BDA0001964567660000064
覆盖周期T表明了当前轨道参数下，交线真近点角需多久时间能变化360°范围。
4.停泊轨道参数设计
由上述分析可以知道，由于摄动的关系，升交点赤经和近地点角会随时间发生变化，升交点赤经的变化又使得交线矢量发生偏转，而交线矢量和近地点矢量的偏转又共同造成了交线地心距差ΔR0和ΔR180随时间发生变化，图1较形象地展示交线地心距差与交线矢量偏转与目标轨道和停泊轨道的近地点矢量偏转的关系。当停泊轨道选择合适的参数时，能使得ΔR0或ΔR180在摄动作用下就能够较快速的接近甚至相交。下面分析停泊轨道在什么条件能最大可能的接近大部分低轨道碎片分析，若是有目标轨道具体范围，也可以依照摄动力影响最大原则设计合适的停泊轨道。
(1)高度选择
LEO区域的空间碎片在轨道高度700～1100km有最大分布，停泊轨道一方面要能覆盖上述轨道高度，另从式(4)(5)可以知道，在相同半长轴情况下，偏心率越大，升交点赤道和近地点角变化越快，当然也需避免轨道高度太低时因大气阻力会造成卫星的轨道衰减过快，最终选择近地点高度400km，远地点高度在1400km的椭圆轨道作为停泊轨道，因此停泊轨道半长轴：
Figure BDA0001964567660000065
停泊轨道偏心率：
Figure BDA0001964567660000066
(2)倾角选择
由式(10)可以看出，
Figure BDA0001964567660000067
与
Figure BDA0001964567660000068
和
Figure BDA0001964567660000069
都有关，
Figure BDA00019645676600000610
项影响因素比较多，与目标轨道参数和停泊轨道参数都有关，而
Figure BDA00019645676600000611
项只跟停泊轨道本身参数有关，因此先从
Figure BDA0001964567660000071
的角度进行分析，将停泊轨道参数：半长轴aP＝7278.14km，偏心率eP＝0.0687，倾角iP范围在0到180度代入式(10)，得近地点变化量
Figure BDA0001964567660000072
随倾角iP的变化情况，具体如图2所示。
当
Figure BDA0001964567660000073
值越大，停泊轨道近地点角变化越快，由图2所示，在iP＝0°或iP＝180°时
Figure BDA0001964567660000074
最大，因此下面专门针对iP＝0°和iP＝180°进行分析。
在iP＝0°和iP＝180°时，对应的交线真近点角fPc可按下式：
Figure BDA0001964567660000075
交线真近点角变化量
Figure BDA0001964567660000076
为：
Figure BDA0001964567660000077
以目标轨道参数偏心率eT＝0，半长轴aT＝7078.14km，目标轨道倾角iT的范围在0～180°，停泊轨道参数：半长轴aP＝7278.14km，偏心率eP＝0.0687，倾角iP＝0°、iP＝180°代入式(13)可得交线真近点角变化率
Figure BDA0001964567660000078
随目标倾角的变化情况，具体如图3所示。
当目标轨道倾角小于等于90°时，停泊轨道倾角iP＝0°，交线真近点角变化率
Figure BDA0001964567660000079
更大，当目标轨道倾角大于90°，停泊轨道倾角iP＝180°，交线真近点角变化率
Figure BDA00019645676600000710
更大。虽然只是eT＝0，aT＝7078.14km的仿真结果，但其它目标轨道也会有类似的结论，当前的碎片绝大部分目标轨道都在99°以下，因此iP＝0°作为停泊轨道设计首选。
下面计算了目标轨道倾角在0和99°时不同轨道高度时交线真近点角变化率
Figure BDA00019645676600000711
和覆盖周期值，由表1、表2、表3和表4可知：当目标轨道倾角iT＝0，目标轨道近地点高度和远地点高度都为400km时，交线真近点角变化率
Figure BDA0001964567660000081
变化最快，为每天西进14.39°，覆盖周期最短为25.02天；当目标轨道倾角iT＝99°，目标轨道近地点高度和远地点高度都为400km时，交线真近点角变化率
Figure BDA0001964567660000082
变化最慢，为每天西进5.08°，覆盖周期最长70.87天。其它目标轨道条件的交线真近点角变化率
Figure BDA0001964567660000083
和覆盖周期都在上述值中间。也就是说，当选择倾角iP＝0，近地点高度400km，远地点高度1400km的停泊轨道时，对于目标轨道倾角在99°以下时，轨道高度在400～1400km，在摄动的作用下，最快只需25.02天，最慢也只需70.87天，fPc就能覆盖一遍360°，从而使得交线地心距差接近0。
表1目标轨道400～1400km，iT＝0，真近点角变化率
Figure BDA0001964567660000084
(°/天)
Figure BDA0001964567660000085
表2目标轨道400～1400km，iT＝99°，真近点角变化率
Figure BDA0001964567660000086
(°/天)
Figure BDA0001964567660000087
表3目标轨道400～1400km，iT＝0，覆盖周期T(天)
Figure BDA0001964567660000091
表4目标轨道400～1400km，iT＝99°，覆盖周期T(天)
Figure BDA0001964567660000092
(二)最优轨控点
从精确外推轨道数据里插值得到目标轨道交线幅角uT＝uTc或uT＝uTc+180°时对应的时刻及停泊轨道和目标轨道数据，将数据代入式(3)(2)对应的交线地心距离差
Figure BDA0001964567660000093
和
Figure BDA0001964567660000094
变化曲线，式(2)只是用二体模型计算，考虑摄动对地心距的影响不超过20km，当目标轨道地心距高度都在[400km,1400km]范围内时，取出
Figure BDA0001964567660000095
和
Figure BDA0001964567660000096
的数据，针对每一个被选出的点，按照如下步骤计算出速度增量最小值min(Δvmin)：
1.可能交会点和可能轨控点
将
Figure BDA0001964567660000097
和
Figure BDA0001964567660000098
对应目标轨道的点为可能交会点，标记为B点，记录的时刻为tB，可能交会点的交线幅角为uTc，轨道参数为(aTB,eTB,iTB,ΩTB,wTB,uTB)，同时取出同时刻停泊轨道数据：交线幅角为
Figure BDA0001964567660000099
轨道参数
Figure BDA00019645676600000910
从可能交会点对应时刻tB往回推1轨，通过插值得到360个可能轨控点，每两个可能轨控点间隔1度，设可能轨控点标记为A点，可能轨控点的交线幅角为uPc，对应的轨道参数为(aPA,ePA,iPA,ΩPA,wPA,uPA)。另外，还需插值得到停泊轨道在交线矢量方向上的点，记为(aPB,ePB,iPB,ΩPB,wPB,uPB)，其中，uPB＝uPc+uTB-uTc，a为半长轴，e为偏心率，i为倾角，Ω为升交点赤经，w为近地点幅角，u为纬度幅角。下标TB表示可能交会点，下标P,t＝tB为停泊轨道在时间tB上对应的点，下标PA表示可能轨控点，下标PB为停泊轨道在交线方向上与目标轨道交会点对应的点。
2.平面投影法后的位置和速度
为避免由于轨道面偏转带来的额外速度增量，在此采用平面投影法，即将目标轨道和停泊轨道都投影到倾角为0，升交点为0的虚拟平面上，然后算出平面投影后的轨道位置和速度，具体步骤如下：
平面投影后可能轨控点的位置矢量为：
Figure BDA0001964567660000101
其中，
Figure BDA0001964567660000102
交线幅角差ΔuPA＝uPA-uPc，真近点fPA＝uPA-wPA。
可能轨控点径向速度：
Figure BDA0001964567660000103
其中，μ＝398600.4418km3/sec2。
可能轨控点飞行方向速度：
Figure BDA0001964567660000104
平面投影后的可能轨控点的速度矢量为：
Figure BDA0001964567660000105
平面投影后可能交会点的位置矢量为：
Figure BDA0001964567660000111
其中，
Figure BDA0001964567660000112
ΔuTB＝uTB-uTc，fTB＝uTB-wTB。
3.虚拟交会点
求解可能轨控点与可能交会点间的最小速度增量时，虽然已经用平面投影法克服法向偏差，但可能轨控点用二体模型外推到可能交会点方向和精确模型外推到交会点方向时两者还是会有径向偏差，这也是由于摄动带来的。而Lagrange转移时间方程只针对二体情况建立两点转移方程，要克服摄动带来的径向偏差，通过建立虚拟交会点后与可能轨控点求解两点最小速度增量，所谓虚拟的交会点就是事先计算出二体模型与精确模型外推后的偏差矢量，然后将此偏差叠加到当前交会点后形成的交会点，具体计算如下：
可能轨控点当前的交线幅角差为ΔuPA，交会点对应的交线幅角差为ΔuB，则两者角度差为Δf0＝ΔuTB-ΔuPA，从可能轨控点外推到可能交会点方向时的径向距离为：
Figure BDA0001964567660000113
可能轨控点A用二体模型外推到交会点方向B时的位置矢量为：
Figure BDA0001964567660000114
其中，
Figure BDA0001964567660000115
计算可能轨控点A到交线幅角差在考虑各种摄动情况下外推到交会点方向时的位置
Figure BDA0001964567660000116
Figure BDA0001964567660000117
其中，
Figure BDA0001964567660000118
fPB＝uPB-wPB。
两者差为：
Figure BDA0001964567660000121
虚拟交会点为：
Figure BDA0001964567660000122
虚拟交会点的模为
Figure BDA0001964567660000123
4.转移角
转移轨道的法向单位矢量：
Figure BDA0001964567660000124
停泊轨道的法向单位矢量：
Figure BDA0001964567660000125
停泊轨道与转移轨道的夹角为：
Figure BDA0001964567660000126
依据速度增量小原则求解转移角Δf以及轨道面的夹角余弦cosε，具体如下：
Figure BDA0001964567660000127
5.求解速度增量最小的转移轨道
在求解最小速度增量最小的转移轨道时，由于要用到Lagrange转移时间方程，以及控后速度、速度增量、速度增量与转移半长轴等知识点，在此先阐述这里理论。
■Lagrange转移时间方程输入参数：
两点弦长c：
Figure BDA0001964567660000128
周长之半s：
Figure BDA0001964567660000129
最小能量半长轴am：
Figure BDA00019645676600001210
■Lagrange转移时间方程及α、β参数：
Lagrange转移时间方程为：
ntr(tB-tA)＝(α-β)-(sinα-sinβ) (31)
式中，ntr为转移轨道角速度，
Figure BDA0001964567660000131
atr为转移轨道半长轴。
记α0和β0为余弦主值，具体表示如下：
Figure BDA0001964567660000132
Figure BDA0001964567660000133
其中，0≤α0≤π，0≤β0≤βm＜π，
Figure BDA0001964567660000134
α和β依据象限判定如下：
Δf≤π且为第1类轨道转移，α＝α0，β＝β0；
Δf≤π且为第2类轨道转移，α＝2π-α0，β＝β0；
Δf＞π且为第1类轨道转移，α＝α0，β＝-β0；
Δf＞π且为第2类轨道转移，α＝α0，β＝-β0。
■当转移半长轴已知，控后速度、飞行路径角求解如下：
●当转移角Δf≠π：
A点控后速度VtrA为：
Figure BDA0001964567660000135
转移轨道的半通径ptr为：
Figure BDA0001964567660000136
etrsinftrA为：
Figure BDA0001964567660000137
转移轨道偏心率为：
Figure BDA0001964567660000138
最终飞行路径角γtrA：
Figure BDA0001964567660000141
●当转移角Δf＝π：
控后速度同(34)，转移轨道的半通径ptr
Figure BDA0001964567660000142
etr sin ftrA为：
Figure BDA0001964567660000143
最终飞行路径角γtrA：
Figure BDA0001964567660000144
■速度增量ΔvA求解：
Figure BDA0001964567660000145
■速度增量ΔvA与转移半长轴atr的导数
Figure BDA0001964567660000146
求解：
当Δf≠π时，
Figure BDA0001964567660000147
当Δf＝π时，
Figure BDA0001964567660000148
在
Figure BDA0001964567660000151
和
Figure BDA0001964567660000152
及转移角Δf确定后，可求解出所有可能轨控点和交会点之间的最小速度增量，然后对最小速度增量统计，统计出速度增量最小的轨控点和对应的轨道，每一个可能轨控点与交会点的求解最小速度增量的步骤如下：
1)计算Lagrange转移时间方程的输入参数：s、c、am。
2)确定转移半长轴区间值，已知可能转移轨道的半长轴范围在[0.5aP 1.5aP],若0.5aP≤am，半长轴区间下限值为ax＝am，反之半长轴计算区间值下限值为ax＝0.5aP，若1.5aP≤am，则当前区间值无解，反之半长轴上限值为as＝1.5aP。
3)当确定转移轨道半长轴区间[ax,as]，按照式(32)～(44)求解
Figure BDA0001964567660000153
和
Figure BDA0001964567660000154
然后根据
Figure BDA0001964567660000155
和
Figure BDA0001964567660000156
的符号特性进行求解轨控点和目标点的最小速度增量转移轨道。
●当第1类轨道转移
Figure BDA0001964567660000157
第2类轨道转移的
Figure BDA0001964567660000158
则最小速度增量发生在atr＝ax时。
●当第1类轨道转移
Figure BDA0001964567660000159
第2类轨道转移的
Figure BDA00019645676600001510
最小速度增量发生在第2类轨道转移：若
Figure BDA00019645676600001511
第2类轨道转移的速度增量还是保持随半长轴的增加继续减小的特性，此时最小速度增量发生在atr＝as时，若
Figure BDA00019645676600001512
最小速度增量发生在
Figure BDA00019645676600001513
用对分法求解。
●当第1类轨道转移
Figure BDA00019645676600001514
第2类轨道转移的
Figure BDA0001964567660000161
则最小速度增量在第1类轨道转移：若
Figure BDA0001964567660000162
第1类轨道转移的速度增量还是保持随半长轴的增加继续减小的特性，此时最小速度增量在atr＝as时，若
Figure BDA0001964567660000163
最小速度增量发生在
Figure BDA0001964567660000164
用对分法求解。
对分法的求解过程为：
i.已知axi、asi，令azi＝(axi+asi)/2，计算
Figure BDA0001964567660000165
其中i初值从0开始迭代，ax0＝ax,as0＝as；
ii.若
Figure BDA0001964567660000166
或|asi-axi|＜eps，停止计算，取atr＝azi，求取最小速度增量值
Figure BDA0001964567660000167
iii.若
Figure BDA0001964567660000168
取ax(i+1)＝axi,as(i+1)＝azi；若
Figure BDA0001964567660000169
取ax(i+1)＝azi,as(i+1)＝asi，转i继续求解。
按照上述过程，360个可能轨控点与对应的可能交会点都能求解出最小速度增量，然后对360个最小速度增量统计，得出速度增量最小的值，记为min(Δvmin)，按照上述计算过程，对于每个可能交会点，都会有对应的最优轨控点和速度增量最小值。
6.最佳交会点、最佳轨控点
当目标轨道在可能交会点上时，停泊轨道当前位置距目标交会点的角
Figure BDA00019645676600001610
由于相位
Figure BDA00019645676600001611
有大有小，对应调相的速度增量有正有负，调相时间有长有短，因此需要考虑
Figure BDA00019645676600001612
和
Figure BDA00019645676600001613
两个方向的求解，当求出
Figure BDA00019645676600001614
和
Figure BDA00019645676600001615
集合中的每一点对应的速度增量为min(Δvmin)及对应的轨控点，无论是那个方向调相，都在对应的轨控点计划两次控制，首次主要功能为调相，第二次主要功能为接近，两个方向的初算出速度增量总和：
Figure BDA0001964567660000171
其中，aP为停泊轨道半长轴，Δt为对应的调相时间，要使得轨控点小，要尽量使得调相时间长。
最终能得到集合
Figure BDA0001964567660000172
然后计算出速度增量总和最小的点：
Figure BDA0001964567660000173
其对应的可能交会点和可能轨控点为最佳交会点和最佳轨控点。
7.迭代计算转移轨道
当选出
Figure BDA0001964567660000174
对应在
Figure BDA0001964567660000175
方向时，则取出对应的速度增量作为第一次控制时的速度增量初值(LVLH坐标系，下同)：
Figure BDA0001964567660000176
若是对应在
Figure BDA0001964567660000177
方向时，则取出对应的速度增量作为第一次控制时的速度增量初值，
Figure BDA0001964567660000178
第二次控制对应的速度增量初值Δv2＝min(Δvmin)-Δv1。
当机动卫星与目标碎片的接近距离要求矢量为d，将此矢量叠加到交会点上作为新的交会点，然后在最佳轨控点代入上述两次速度增量，得到过渡闭合轨道参数和转移轨道参数，然后按照下列步骤进行迭代计算最终的最优转移轨道。
1)取出过渡闭合轨道半长轴
Figure BDA0001964567660000179
转移轨道半长轴
Figure BDA00019645676600001710
依据下式进行计算：
Figure BDA00019645676600001711
Figure BDA00019645676600001712
其中，
Figure BDA0001964567660000181
Figure BDA0001964567660000182
是上一次迭代后到对应交会时刻机动卫星在转移轨道的相位。
2)然后更新第一次速度增量：
Figure BDA0001964567660000183
将更新后的Δv1代入模型继续外推，Δv2不变，取出第二次控后的数据到交会时刻的数据，依照平面投影法计算出以交线为起点的坐标系下的位置和数据，然后求解目标星在交会时刻时与机动卫星的距离，其值小于1m时，停止迭代，否则继续按照单脉冲轨道设计的方法，先求出虚拟交会点，然后求出第二次轨控结束时刻的位置和速度，依据此两点求解最小速度增量
Figure BDA0001964567660000184
虚拟交会点和第二次轨控结束时刻的点依然采用平面投影后的位置和速度，最后更新第二次速度增量：
Figure BDA0001964567660000185
将更新后的速度增量Δv2代入模型继续精确外推，继续1)循环。
以下是基于本实施例设计方法的一个仿真实验及其结果。
目标：机动卫星与CZ-2C火箭碎片最近的距离100m，无方向要求。
(一)卫星参数
机动卫星参数：
(1)卫星初始质量：100kg
(2)迎风面积：0.3m2
(3)辐射面积：0.3m2
(4)推进参数：
推力：20N
单次最长工作时长：1200s
单次提供的最大速度增量：约240m/s
初始燃料：31.5kg
比冲：2156Ns/kg
(5)最小轨控间隔：7h
目标碎片参数：
质量：1000kg
迎风面积：2m2
辐射面积：2m2
目标是2018年发射的国际编号为43172的CZ-2C火箭碎片，初始轨道起始时间为：29Jan 2018 13:41:06.000UTCG，其他轨道参数如下：
表5初始轨道参数(瞬根数)
参数 机动卫星 CZ-2C 差值
半长轴(km) 7278.137 6983.084 -295.053
偏心率 0.068699 0.001974 -0.066725
倾角(°) 0 35.048 35.048
升交点赤经(°) 0 158.457 158.457
近地点幅角(°) 50 351.817 301.817
纬度幅角(°) 180 0.202 -179.798
真近点角(°) 0 8.385 8.385
平近点角(°) 0 8.352 8.352
(二)仿真工具与模型
整个仿真工具及参数设置、仿真工况如下：
(1)仿真工具
STK 10：Astrogator模型
Matlab2010b
(2)模型参数
引力模型采用WGS84_EGM96，参数设置：Degree:70，Order:70。
大气模型采用MSISE 1990，参数按照stk默认设置。
光压模型采用Spherical SRP，参数按照stk默认设置。
第三体引力考虑太阳和月亮。
机动卫星与目标卫星的大气阻力系数Cd为2.2。
机动卫星与目标卫星的辐射系数Cr为1。
(三)标称仿真结果
覆盖周期29.32天，一个覆盖周期交线地心距差随时间的变化情况见图4。
选取小于20km的交线地心距差，求取对应的最小速度增量，以及调相速度增量，求解两个方向的速度总和
Figure BDA0001964567660000191
然后求出
Figure BDA0001964567660000192
具体如下表格：
表6相关参数计算结果
Figure BDA0001964567660000193
Figure BDA0001964567660000201
上述表格中标灰的分别是速度增量最小和次小的点，这两个点的速度增量接近，但第一个点的时间确比第二点的时间小14天，因此选第一个点继续后续的仿真。然后通过matlab和stk的联合仿真，可以得到：经过2次地面，控制耗时13.5天，机动卫星在12Feb2018 00:57:09.678UTCG与CZ-2C碎片距离99.2m的地方飞过，总共需速度增量1.039m/s，整体交会图见图5，其中标号，1表示目标轨道、2表示停泊轨道、3表示闭合轨道、4表示交会轨道、5表示交会点，具体控制参数和仿真结果如下：
表7控制参数(第3次之后是星上自主控制)
Figure BDA0001964567660000202
从图6可以看出，机动卫星第13次接近时，以距离99.57m飞越CZ-2C，与目标距离100m的误差为0.43m。
本实施例的设计方法，依据近地轨道的摄动特性，设计出相应的椭圆轨道，首先利用摄动，实现停泊轨道与目标轨道两者在轨道面交线方向上距离的接近，然后采用卫星控制，实现机动卫星与目标碎片的接近。通过摄动主控，卫星辅控的方式，只需很少的燃料就能接近碎片，在相同的燃料下能完成更多的接近动作。上述设计方法适用于接近碎片，以及其它对时间要求不高的航天任务。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (9)

1.一种用于空间碎片接近的轨道设计方法，其特征在于，包括：
设定机动卫星的停泊轨道；
将所述停泊轨道和目标轨道的交线地心距差足够小的交线点取出作为可能交会点，根据所述可能交会点取出对应的可能轨控点；
在无时间的约束下，求解每个所述可能交会点和与其相应的所述可能轨控点之间的最小速度增量，从而每个所述可能交会点均可得到与其对应的速度增量最小值和最优轨控点；
结合相位调整对速度增量的影响，优选出最佳交会点和最佳轨控点；
最后将所述机动卫星与目标碎片在接近时的距离矢量要求叠加到所述最佳交会点上，迭代计算最终所需速度增量，得到最优转移轨道。
2.如权利要求1所述的设计方法，其特征在于，所述停泊轨道被设置为椭圆轨道。
3.如权利要求2所述的设计方法，其特征在于，所述停泊轨道的参数设置包括：停泊轨道倾角0°、近地点高度400km、远地点高度1400km。
4.如权利要求1所述的设计方法，其特征在于，所述交线地心距差是停泊轨道交点径向地心距与目标轨道交点径向地心距的差异。
5.如权利要求4所述的设计方法，其特征在于，将所述交线地心距差小于20km的交线点取出作为所述可能交会点。
6.如权利要求1所述的设计方法，其特征在于，取出距所述可能交会点最近的一轨数据为所述可能轨控点。
7.如权利要求6所述的设计方法，其特征在于，从每个所述可能交会点对应时刻回推一轨，通过插值得到360个所述可能轨控点，每两个所述可能轨控点间隔1度。
8.如权利要求1所述的设计方法，其特征在于，在求解每个所述可能交会点和与其相应的所述可能轨控点之间的最小速度增量时，通过平面投影法克服所述停泊轨道和所述目标轨道的法向偏差。
9.如权利要求1所述的设计方法，其特征在于，在求解每个所述可能交会点和与其相应的所述可能轨控点之间的最小速度增量时，所述可能交会点采用相应的虚拟交会点，所述虚拟交会点是事先计算出二体模型与精确模型外推后的偏差矢量，然后将所述偏差矢量叠加到当前交会点后形成的交会点。
CN201910096002.8A 2019-01-31 2019-01-31 一种用于空间碎片接近的轨道设计方法 Active CN109774974B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910096002.8A CN109774974B (zh) 2019-01-31 2019-01-31 一种用于空间碎片接近的轨道设计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910096002.8A CN109774974B (zh) 2019-01-31 2019-01-31 一种用于空间碎片接近的轨道设计方法

Publications (2)

Publication Number Publication Date
CN109774974A CN109774974A (zh) 2019-05-21
CN109774974B true CN109774974B (zh) 2020-08-14

Family

ID=66503956

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910096002.8A Active CN109774974B (zh) 2019-01-31 2019-01-31 一种用于空间碎片接近的轨道设计方法

Country Status (1)

Country Link
CN (1) CN109774974B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111102982B (zh) * 2019-12-11 2021-09-24 上海卫星工程研究所 高轨目标的抵近方法
CN111619828B (zh) * 2020-04-20 2021-12-07 中国卫通集团股份有限公司 一种同步轨道卫星离轨的方法及装置
CN112000121B (zh) * 2020-07-14 2022-07-12 哈尔滨工业大学 一种多服务飞行器空间在轨服务燃料最优轨道的设计方法
CN112623277B (zh) * 2020-12-25 2022-04-12 中国空间技术研究院 一种对空间异面圆轨道目标快速抵达的变轨方法
CN113848567B (zh) * 2021-08-26 2023-05-30 深圳市魔方卫星科技有限公司 一种sar卫星面内最优轨控确定方法、装置及相关设备
CN114383619B (zh) * 2021-12-07 2023-09-05 上海航天控制技术研究所 一种高精度轨道计算方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101214859A (zh) * 2007-12-26 2008-07-09 北京控制工程研究所 一种变轨期间自主故障检测恢复控制的方法
CN101916114A (zh) * 2010-04-14 2010-12-15 清华大学 一种为卫星星座提供在轨服务的服务轨道设计方法
CN102999616A (zh) * 2012-11-29 2013-03-27 北京理工大学 一种基于轨道根数的星际飞行发射机会搜索方法
CN103594749A (zh) * 2013-10-30 2014-02-19 中国运载火箭技术研究院 基于无线能源传输的空间飞行器充电方法
CN103970142A (zh) * 2013-02-01 2014-08-06 上海新跃仪表厂 在轨拖曳的组合体航天器姿轨复合控制方法
US10054449B2 (en) * 2014-02-28 2018-08-21 Thales Method of following a transfer orbit or a phase of orbital placement of a space vehicle, in particular an electric propulsion vehicle, and apparatus for the implementation of such a method

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101214859A (zh) * 2007-12-26 2008-07-09 北京控制工程研究所 一种变轨期间自主故障检测恢复控制的方法
CN101916114A (zh) * 2010-04-14 2010-12-15 清华大学 一种为卫星星座提供在轨服务的服务轨道设计方法
CN102999616A (zh) * 2012-11-29 2013-03-27 北京理工大学 一种基于轨道根数的星际飞行发射机会搜索方法
CN103970142A (zh) * 2013-02-01 2014-08-06 上海新跃仪表厂 在轨拖曳的组合体航天器姿轨复合控制方法
CN103594749A (zh) * 2013-10-30 2014-02-19 中国运载火箭技术研究院 基于无线能源传输的空间飞行器充电方法
US10054449B2 (en) * 2014-02-28 2018-08-21 Thales Method of following a transfer orbit or a phase of orbital placement of a space vehicle, in particular an electric propulsion vehicle, and apparatus for the implementation of such a method

Also Published As

Publication number Publication date
CN109774974A (zh) 2019-05-21

Similar Documents

Publication Publication Date Title
CN109774974B (zh) 一种用于空间碎片接近的轨道设计方法
CN109625323B (zh) 一种卫星化学推进变轨方法及系统
EP1654159B1 (en) Apparatus for a geosynchronous life extension spacecraft
Roger et al. Safety constrained free-flyer path planning at the international space station
US9428285B2 (en) System and method for managing momentum accumulation
CN105511493A (zh) 一种基于火星大气辅助的低轨星座部署方法
Baranov et al. Optimal transfer schemes between space debris objects in geostationary orbit
Wang et al. Assembled kinetic impactor for deflecting asteroids by combining the spacecraft with the launch vehicle upper stage
Sah et al. Design development of debris chaser small satellite with robotic manipulators for debris removal
Leomanni et al. All-electric spacecraft precision pointing using model predictive control
Folta et al. Design and implementation of satellite formations and constellations
Malyshev et al. Orbital corrections of space vehicles while performing dynamic operations
Lin et al. Entire flight trajectory design for temporary reconnaissance mission
Xie et al. Guidance, navigation, and control for spacecraft rendezvous and docking: theory and methods
Li et al. Analytical design methods for determining Moon-to-Earth trajectories
Ravikumar et al. Autonomous terminal maneuver of spacecrafts for rendezvous using model predictive control
Ortega Fuzzy logic techniques for rendezvous and docking of two geostationary satellites
Sostaric et al. Lunar ascent and rendezvous trajectory design
Schäff et al. Advanced Electric Orbit-Raising Optimization for Operational Purpose
Biesbroek The e. Deorbit CDF study
Antreasian et al. Preliminary planning for NEAR's low-altitude operations at 433 Eros
Moessner et al. CAT Differential Drag Implementation and Lessons Learned
Rana et al. Study into In-orbit Servicing of the Rosetta Mission
Oki et al. Satellite System Design with The Blocking Effect: Application to Active Debris Removal Mission
Colagrossi et al. Enhancing Technologies and Operations for Service Transportation in Cislunar Environment

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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20230906

Address after: 201306 building C, No. 888, Huanhu West 2nd Road, Lingang New District, China (Shanghai) pilot Free Trade Zone, Pudong New Area, Shanghai

Patentee after: Shanghai Zhongkechen New Satellite Technology Co.,Ltd.

Address before: No. 4 Building, 99 Haike Road, Pudong New Area, Shanghai, 201203

Patentee before: SHANGHAI ENGINEERING CENTER FOR MICROSATELLITES