一种航天器多约束飞越最优速度的估计方法

1.本发明涉及航天器飞越状态的估计技术领域，更具体地说，特别涉及一种航天器多约束飞越最优速度的估计方法。


背景技术：

2.随着星链卫星等低轨巨型星座的建造，近地空间的在轨目标急剧增多，遮天蔽日的卫星群对未来空间态势感知能力提出了严峻的挑战。由于天基光学成像技术的发展，基于装有光学传感器的卫星的空间态势感知体系得到了进一步的重视。利用多个航天器进行异面飞越，可在较短时间内以较小的燃料代价抵近飞越并拍摄多颗目标卫星，从而快速获得星座、编队等集群卫星的图像资料，为后续的任务决策提供支撑。
3.在研究一对多、多对多飞越巡察序列问题时，需要高效的单对单飞越成本估计进行支撑。在交会估计问题中，由于始末两端的位置速度已知，可以根据轨道根数差来估计相应时间下的转移成本。然而，在飞越任务中由于末端速度未知，需要估计出满足终端速度约束的最优飞越速度，使得整个飞越转移成本最小，因此在飞越终端速度不能提前确定情况下的估计飞越任务成本是目前亟待解决的技术问题。


技术实现要素：

4.本发明的目的在于提供一种航天器多约束飞越最优速度的估计方法，以克服现有技术所存在的缺陷。
5.为了达到上述目的，本发明采用的技术方案如下：
6.一种航天器多约束飞越最优速度的估计方法，包括以下步骤：
7.s1、建立飞越巡察的终端不等式约束表达式；
8.s2、通过航天器与目标的轨道异面度计算飞越最优时刻和飞越轨道面方位；
9.s3、根据飞越轨道面方位给出飞越速度的初值，基于库恩-塔克条件建立固定时间下的多约束飞越最优速度估计。
10.进一步地，所述步骤s1具体为：
11.设定在地心惯性坐标系下航天器起始时刻t0和飞越时刻t的位置速度为[r
c0
,v
c0
]和[rc,vc]，以及目标在起始时刻t0和飞越时刻t的位置速度为[r
t0
,v
t0
]和[r
t
,v
t
]，设定飞越任务为位置重合、速度存在偏差的无偏飞越；设定飞越前的设定时间段航天器和目标的运动均是匀速直线的，由此可以推得航天器和目标相对位置矢量δr＝r
t
–
rc与相对速度矢量δv＝v
t
–
vc的叉乘为0，并用δv来代替位置重合前的视线方向；设定视线方向与太阳位置矢量rs方向的夹角下限为α
smin
，得到关于视线阳光角的约束表达式为：
[0012]
arccos((rs/rs)
·
(δv/δv))-(π-α
smin
)≤0
ꢀꢀ
(36)
[0013]
设定航天器飞越时的相对速度大小不能超过上限δv
max
：
[0014]
δv-δv
max
≤0
ꢀꢀ
(37)
[0015]
设定航天器飞越时的近地点高度不能低于下限r
pmin
：
[0016]rpmin-ac(1-ec)≤0
ꢀꢀ
(38)
[0017]
式中，ac和ec为航天器在飞越时刻t的半长轴和偏心率。
[0018]
进一步地，所述步骤s2具体包括以下步骤：
[0019]
s20、设定航天器在起始时刻t0和飞越时刻t的轨道根数为[a
c0
,e
c0
,i
c0
,ω
c0
,ω
c0
,m
c0
]和[ac,ec,ic,ωc,ωc,mc]，以及目标在起始时刻t0和飞越时刻t的轨道根数为[a
t0
,e
t0
,i
t0
,ω
t0
,ω
t0
,m
t0
]和[a
t
,e
t
,i
t
,ω
t
,ω
t
,m
t
]，a,e,i,ω,ω,m分别为半长轴、偏心率、倾角、升交点赤径、近拱点纬度幅角和平近点角，下标c0和c为航天器在t0和t时刻的轨道根数，下标t0和t为目标在t0和t时刻的轨道根数，轨道转移时长δt＝t
–
t0，根据球面三角形公式计算出两个初始轨道面的异面度φ以及相对速度约束所能容许的最大异面度φ
max
为：
[0020][0021]
式中，μ为中心天体的万有引力常数；
[0022]
s21、当φ≤φ
max
，航天器允许只通过共面机动来飞越目标，此时最优飞越时刻δt
opt
的表达式为：
[0023][0024]
式中，u
t
为航天器和目标的初始轨道面交点在目标轨道上的相位，u
t0
为目标的初始相位，k为非负整数，u
t
的表达式为：
[0025][0026]
式中，u
t
有两个解，只取不在地影一侧的解，和为目标在仅考虑j2长期项摄动下的轨道根数进动率，包括航天器在内的轨道根数进动率表达式为：
[0027][0028]
式中，j2为摄动项常数，re为地球平均赤道半径；
[0029]
s22、根据φ与φ
max
，以及δt与δt
opt
的关系，确定飞越轨道面的单位法向量矢量he，记时间偏差的阈值为δt
threshold
，当φ≤φ
max
并且|δt-δt
opt
|≤δt
threshold
，认为航天器可以进行共面机动飞越，则根据球面三角形中的角角边(i
c0-i
t0-u
t
)来求解升交点赤径差δω为：
[0030]
[0031]
当φ≤φ
max
并且|δt-δt
opt
|》δt
threshold
，或者φ》φ
max
时，航天器需要异面机动变轨，根据球面三角形中的角边角(i
t0-u
t-φ
max
)来求解倾角ic和升交点赤径差δω为：
[0032][0033]
根据δω＝ω
t0
–
ωc，由此得到单位法向量he为
[0034][0035]
进一步地，所述步骤s3具体包括以下步骤：
[0036]
s30、建立以目标为原点o，r
t
方向为x轴，he方向为z轴的lvlh相对坐标系，设定飞越速度的大小和三个方向分量为vc,γ
x
,γy,γz，速度矢量vc与he正交，飞越速度的大小vc在xoy平面内，设定飞越速度的大小vc与y轴正向的夹角为θ∈[-θ
max
,θ
max
]，θ
max
为角度的边界值，通过坐标系变化，可用θ描述速度的三个方向分量为：
[0037][0038]
式中，为lvlh坐标系到地心惯性坐标系的旋转矩阵；
[0039]
速度的大小根据位置速度与轨道根数的置换公式推得关于偏心率ec的表达式：
[0040][0041]
式中，e
cmax
是飞越轨道偏心率的容许上限；
[0042]
通过离散vc和θ在各自上下界中的取值，得猜测值为：
[0043][0044]
式中，n1和n2为离散点的数量，α和β均为整数，并将θ
guess
带入式(11)中得方向的猜测初值；
[0045]
s31、设定飞越轨迹通过单脉冲进行估计，则其引起的轨道根数偏差为：
[0046]
[0047]
则单次脉冲的速度增量大小δv为：
[0048][0049]
式中，为航天器初始时刻的平均轨道速度；
[0050]
以速度增量的平方为最优问题的目标函数，表达式为：
[0051][0052]
由于速度的三个单位方向的模值为1，引入等式约束为：
[0053][0054]
建立拉格朗日方程，表达式为：
[0055]
l＝f λ
t
g κ
t
σ
ꢀꢀ
(53)
[0056]
式中，λ为等式约束乘子，κ为不等式约束乘子，根据库恩-塔克条件，对四个未知数求偏导，可得到一个8维的等式方程组为：
[0057][0058]
式中，涉及到ac,ec,ic,ωc,ωc对vc,γ
x
,γy,γz的偏导，需要将这些轨道根数用位置速度来表示，记hc＝[h
cx
,h
cy
,h
cz
]
t
，ec＝[e
cx
,e
cy
,e
cz
]
t
，表达式为：
[0059][0060]
其中，hc为飞越轨道的角动量，记γ＝[γ
x
,γy,γz]
t
，可推得偏导数为：
[0061][0062][0063][0064]
[0065][0066][0067][0068]
记i
x
＝[1,0,0]
t
,iy＝[0,1,0]
t
,iz＝[0,0,1]
t
，角动量矢量hc及其模值hc关于vc和γ的偏导为：
[0069][0070][0071]
相对速度矢量δv及其模值δv关于vc和γ的偏导为：
[0072][0073]
偏心率矢量ec关于vc和γ的偏导为：
[0074][0075]
再根据步骤s3中给出的关于飞越速度的初值[vc,γ
x
,γy,γz]0，以及等式约束和不等式约束乘子的初值[λ,κ1,κ2,κ3]0＝0，基于powell的混合数值迭代器进行迭代可得到最终的飞越速度解；
[0076]
s32、得到最终的飞越速度解之后，会产生一个相位差
[0077][0078]
其中，u
cf
为航天器无机动外推至飞越时刻的相位，和为航天器有机动相对于无机动时的进动率偏差，对相对平均角速度进行一阶近似，可得到：
[0079][0080]
利用半长轴偏差进行相位修正，使得δu δncδt＝0，修正后的半长轴表达式为：
[0081][0082]
其中，n∈[-1,1]表示圈数，修正后的飞越速度大小为：
[0083][0084]
最终得到的飞越速度为
[0085]
s33、飞越成本估计转换成交会成本估计，并根据交会估计算法进行估计；
[0086]
s34、重复步骤s30-s33，统计每次迭代生成的速度增量δv，选择最小的值记为δv
best
，对应的飞越速度解为[vc,γ
x
,γy,γz]
best
。
[0087]
与现有技术相比，本发明的优点在于：本发明提供的一种航天器多约束飞越最优速度的估计方法，将末端飞越速度的大小和三个方向分量设为4个未知数，设计优化问题的目标函数为单脉冲下的速度增量，引入包括飞越巡察阳光角、飞越相对速度以及飞越时的轨道近地点高度在内的3个不等式约束，利用库恩-塔克条件列出了一个8维的等式方程组，通过方程组的求解获得飞越最优速度估计，将飞越估计变为交会估计，从而利用已知的轨道根数差来估计转移成本。
附图说明
[0088]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0089]
图1是本发明中航天器飞越时刻的视线阳光角约束示意。
[0090]
图2是本发明航天器多约束飞越最优速度的估计方法流程图。
具体实施方式
[0091]
下面结合附图对本发明的优选实施例进行详细阐述，以使本发明的优点和特征能更易于被本领域技术人员理解，从而对本发明的保护范围做出更为清楚明确的界定。
[0092]
参阅图1和图2所示，本实施例公开了一种航天器多约束飞越最优速度的估计方法，包括以下步骤：
[0093]
步骤s1、建立飞越巡察的终端不等式约束表达式。
[0094]
具体的，设定在地心惯性坐标系下航天器起始时刻t0和飞越时刻t的位置速度为[r
c0
,v
c0
]和[rc,vc]，以及目标在起始时刻t0和飞越时刻t的位置速度为[r
t0
,v
t0
]和[r
t
,v
t
]，设定飞越任务为位置重合、速度存在偏差的无偏飞越；设定飞越前的设定时间段航天器和目标的运动均是匀速直线的，由此可以推得航天器和目标相对位置矢量δr＝r
t
–
rc与相对速度矢量δv＝v
t
–
vc的叉乘为0，并用δv来代替位置重合前的视线方向；设定视线方向与太阳位置矢量rs方向的夹角下限为α
smin
，得到关于视线阳光角的约束表达式为：
[0095]
arccos((rs/rs)
·
(δv/δv))-(π-α
smin
)≤0
ꢀꢀ
(71)
[0096]
设定航天器飞越时的相对速度大小不能超过上限δv
max
：
[0097]
δv-δv
max
≤0
ꢀꢀ
(72)
[0098]
设定航天器飞越时的近地点高度不能低于下限r
pmin
：
[0099]rpmin-ac(1-ec)≤0
ꢀꢀ
(73)
[0100]
式中，ac和ec为航天器在飞越时刻t的半长轴和偏心率。
[0101]
步骤s2、通过航天器与目标的轨道异面度计算飞越最优时刻和飞越轨道面方位。
[0102]
具体的，步骤s2包括以下步骤：
[0103]
步骤s20、设定航天器在起始时刻t0和飞越时刻t的轨道根数为[a
c0
,e
c0
,i
c0
,ω
c0
,ω
c0
,m
c0
]和[ac,ec,ic,ωc,ωc,mc]，以及目标在起始时刻t0和飞越时刻t的轨道根数为[a
t0
,e
t0
,i
t0
,ω
t0
,ω
t0
,m
t0
]和[a
t
,e
t
,i
t
,ω
t
,ω
t
,m
t
]，a,e,i,ω,ω,m分别为半长轴、偏心率、倾角、升交点赤径、近拱点纬度幅角和平近点角，下标c0和c为航天器在t0和t时刻的轨道根数，下标t0和t为目标在t0和t时刻的轨道根数，轨道转移时长δt＝t
–
t0，根据球面三角形公式计算出两个初始轨道面的异面度φ以及相对速度约束所能容许的最大异面度φ
max
为：
[0104][0105]
式中，μ为中心天体的万有引力常数。
[0106]
步骤s21、当φ≤φ
max
，航天器允许只通过共面机动来飞越目标，此时最优飞越时刻δt
opt
的表达式为：
[0107][0108]
式中，u
t
为航天器和目标的初始轨道面交点在目标轨道上的相位，u
t0
为目标的初始相位，k为非负整数，u
t
的表达式为：
[0109][0110]
式中，u
t
有两个解，只取不在地影一侧的解，和为目标在仅考虑j2长期项摄动下的轨道根数进动率，包括航天器在内的轨道根数进动率表达式为：
[0111][0112]
式中，j2为摄动项常数，re为地球平均赤道半径。
[0113]
步骤s22、根据φ与φ
max
，以及δt与δt
opt
的关系，确定飞越轨道面的单位法向量矢量he，记时间偏差的阈值为δt
threshold
，当φ≤φ
max
并且|δt-δt
opt
|≤δt
threshold
，认为航天器可以进行共面机动飞越，则根据球面三角形中的角角边(i
c0-i
t0-u
t
)来求解升交点赤径差δω为：
[0114][0115]
当φ≤φ
max
并且|δt-δt
opt
|》δt
threshold
，或者φ》φ
max
时，航天器需要异面机动变轨，根据球面三角形中的角边角(i
t0-u
t-φ
max
)来求解倾角ic和升交点赤径差δω为：
[0116][0117]
根据δω＝ω
t0
–
ωc，由此得到单位法向量he为
[0118][0119]
步骤s3、根据飞越轨道面方位给出飞越速度的初值，基于库恩-塔克条件建立固定时间下的多约束飞越最优速度估计。
[0120]
具体的，步骤s3具体包括以下步骤：
[0121]
步骤s30、建立以目标为原点o，r
t
方向为x轴，he方向为z轴的lvlh相对坐标系，设定
飞越速度的大小和三个方向分量为vc,γ
x
,γy,γz，速度矢量vc与he正交，飞越速度的大小vc在xoy平面内，设定飞越速度的大小vc与y轴正向的夹角为θ∈[-θ
max
,θ
max
]，θ
max
为角度的边界值，通过坐标系变化，可用θ描述速度的三个方向分量为：
[0122][0123]
式中，为lvlh坐标系到地心惯性坐标系的旋转矩阵；
[0124]
速度的大小根据位置速度与轨道根数的置换公式推得关于偏心率ec的表达式：
[0125][0126]
式中，e
cmax
是飞越轨道偏心率的容许上限；
[0127]
通过离散vc和θ在各自上下界中的取值，得猜测值为：
[0128][0129]
式中，n1和n2为离散点的数量，α和β均为整数，并将θ
guess
带入式(11)中得方向的猜测初值；
[0130]
步骤s31、设定飞越轨迹通过单脉冲进行估计，则其引起的轨道根数偏差为：
[0131][0132]
则单次脉冲的速度增量大小δv为：
[0133][0134]
式中，δva=0.5v
c0
δa/a
c0
，δvi=v
c0
δi，δv
ω
=v
c0
δωsini
c0
，v
c0
为航天器初始时刻的平均轨道速度；
[0135]
以速度增量的平方为最优问题的目标函数，表达式为：
[0136][0137]
由于速度的三个单位方向的模值为1，引入等式约束为：
[0138][0139]
建立拉格朗日方程，表达式为：
[0140]
l＝f λ
t
g κ
t
σ
ꢀꢀ
(88)
[0141]
式中，λ为等式约束乘子，κ为不等式约束乘子，根据库恩-塔克条件，对四个未知数求偏导，可得到一个8维的等式方程组为：
[0142][0143]
式中，涉及到ac,ec,ic,ωc,ωc对vc,γ
x
,γy,γz的偏导，需要将这些轨道根数用位置速度来表示，记hc＝[h
cx
,h
cy
,h
cz
]
t
，ec＝[e
cx
,e
cy
,e
cz
]
t
，表达式为：
[0144][0145]
其中，hc为飞越轨道的角动量，记γ＝[γ
x
,γy,γz]
t
，可推得偏导数为：
[0146][0147][0148][0149]
[0150][0151][0152][0153]
记i
x
＝[1,0,0]
t
,iy＝[0,1,0]
t
,iz＝[0,0,1]
t
，角动量矢量hc及其模值hc关于vc和γ的偏导为：
[0154][0155][0156]
相对速度矢量δv及其模值δv关于vc和γ的偏导为：
[0157][0158]
偏心率矢量ec关于vc和γ的偏导为：
[0159][0160]
再根据步骤s3中给出的关于飞越速度的初值[vc,γ
x
,γy,γz]0，以及等式约束和不等式约束乘子的初值[λ,κ1,κ2,κ3]0＝0，基于powell的混合数值迭代器minpack-1【more j.j.,garbow b.s.,and hillstrom k.e.userguide for minpack-1.argonne national lab.rept.anl-80-74,illinois,1980.】进行迭代可得到最终的飞越速度解。
[0161]
步骤s32、得到最终的飞越速度解之后，会产生一个相位差
[0162][0163]
其中，u
cf
为航天器无机动外推至飞越时刻的相位，和为航天器有机动相对于无机动时的进动率偏差，对相对平均角速度进行一阶近似，可得到：
[0164][0165]
利用半长轴偏差进行相位修正，使得δu δncδt＝0，修正后的半长轴表达式为：
[0166][0167]
其中，n∈[-1,1]表示圈数，修正后的飞越速度大小为：
[0168][0169]
最终得到的飞越速度为
[0170]
步骤s33、当得到终端速度后，即飞越成本估计转换成交会成本估计，并根据交会估计算法进行估计，交会估计算法可以参考现有成熟的摄动交会估计算法【罗亚中,黄岸毅,李恒年,伍升钢,张进,杨震.一种j2摄动下长时间轨道交会最优速度增量快速估计方法[p].湖南省：cn110789739b,2020-12-11.】。
[0171]
步骤s34、重复步骤s30-s33，统计每次迭代生成的速度增量δv，选择最小的值记为δv
best
，对应的飞越速度解为[vc,γ
x
,γy,γz]
best
。如流程图2所示。
[0172]
下面通过实施例对本发明作进一步描述。
[0173]
1)建立飞越巡察的终端不等式约束表达式；
[0174]
给定α
smin
＝30deg，给定δv
max
＝2km/s，给定r
pmin
＝6600km。
[0175]
2)通过航天器与目标的轨道异面度计算飞越最优时刻和飞越轨道面方位。
[0176]
2.1)记初始历元为t0＝59274mjd，航天器和目标在初始时刻的轨道根数见表1：
[0177]
表1航天器和目标在t0＝59274mjd的轨道根数
[0178][0179]
根据球面三角形公式计算出两个初始轨道面的异面度φ＝2.349deg，最大异面度φ
max
＝15.298deg。
[0180]
2.2)设飞越圈数k＝11，此时最优飞越时刻δt
opt
＝68731.86s，给定飞越转移时间δt＝68731.86s。
[0181]
2.3)根据φ与φ
max
，以及δt与δt
opt
的关系，确定飞越轨道面的单位法向量矢量he＝[0.12908,-0.52305,0.84247]
t
。
[0182]
3)根据飞越轨道面方位给出飞越速度的初值，基于库恩-塔克条件建立固定时间下的多约束飞越最优速度估计。
[0183]
3.1)给定θ
max
＝1.5deg，e
cmax
＝0.02，可得v
cmax
＝7.57km/s，v
cmin
＝7.42km/s。设θ的取值步长为0.5deg，则n2＝7；设vc的取值步长为0.01km/s，则n1＝15。α和β均取0，可得第一个飞越速度的初值为[vc,γ
x
,γy,γz]0＝[7.425,0.464,0.783,0.415]。
[0184]
3.2)根据建立的8维等数方程组并利用minpack-1求解器对速度初值进行迭代，得到飞越速度的迭代值为[vc,γ
x
,γy,γz]＝[7.500,0.436,0.877,0.204]。
[0185]
3.3)在该速度下产生的相位差为96.085
°
，经过修成得到新的飞越速度小大为vc＝7.446km/s。
[0186]
3.4)得到飞越速度后，根据新的末端飞越状态与航天器初始状态的轨道根数差来估计转移成本，得到相应的速度增量δv＝215.6m/s。
[0187]
3.5)循环取值105次，依次执行步骤3.1)至3.4)，取δv最小的解为最优的飞越速度[vc,γ
x
,γy,γz]
best
和最优速度增量δv
best
。为了对比算法的精度和效率，采用差分进化算法对同样算例进行3脉冲飞越轨迹优化，最终优化与估计的结果与相对误差见表2：
[0188]
表2飞越估计与差分进化算法优化结果对比
[0189][0190]
表中，优化的飞越速度4个变量的相对误差较小，最大不超过5％，对应的速度增量误差不超过8％，因此估计算法的精度是可观的。
[0191]
优化时间上，3脉冲飞越轨迹有9个设计变量，设差分进化算法的种群为108，进化代数为200，并行核数为6，运行20次平均需要15s；而本发明提出的算法，每次循环只需要约0.001s，运行105次大约需要0.1s，比差分进化算法少两个数量级。
[0192]
本发明将末端飞越速度的大小和三个方向分量设为4个未知数，设计优化问题的目标函数为单脉冲下的速度增量，引入包括飞越巡察阳光角、飞越相对速度以及飞越时的
轨道近地点高度在内的3个不等式约束，利用库恩-塔克条件列出了一个8维的等式方程组，通过方程组的求解获得飞越最优速度估计，将飞越估计变为交会估计，从而利用已知的轨道根数差来估计转移成本。
[0193]
虽然结合附图描述了本发明的实施方式，但是专利所有者可以在所附权利要求的范围之内做出各种变形或修改，只要不超过本发明的权利要求所描述的保护范围，都应当在本发明的保护范围之内。