本发明公开的一种J2摄动下长时间轨道交会最优速度增量快速估计方法。首先计算航天器无机动时与目标在交会时刻的轨道根数差；然后计算仅消除升交点赤经差时的最优控制量；将其作为初值，代入同时消除半长轴、倾角和升交点赤经差的控制量优化模型，解得新的最优控制量；然后再考虑根据轨道面内相位差对控制量进行修正；最后再根据偏心率矢量差对控制量进行修正；最终得到估计的最优交会速度增量。本发明公开的方法解决了无需精确求解时快速估计最优交会速度增量的问题，估计精度较高，适用于有初始与目标轨道的升交点赤经差较大，其它轨道根数偏差较小的情形。
一种J2摄动下长时间轨道交会最优速度增量快速估计方法
技术领域
[0001] 本发明属于航天导航控制技术领域，具体涉及一种J2摄动下长时间轨道交会最优速度增量快速估计方法。
背景技术
[0002] 在近地轨道上，J2项是最主要的摄动项。在进行交会轨道优化设计时，必须考虑J2项的影响，但相应的会使计算量比二体动力学模型增加。
[0003] 现有的轨道交会最优速度增量求解方法多为数值优化算法，计算单次多脉冲轨道转移至少需要几分钟甚至几个小时。在某些轨道设计应用中，可能并不关心具体每个脉冲的大小和方向，只需要获知总的速度增量大小，此时就需要有快速准确的估计方法。现有估计方法多数集中在轨道面内转移速度增量的估计，对于轨道面外交会问题缺少相应研究。
发明内容
[0004] 针对现有技术存在的问题，本发明的目的是提供一种J2摄动下长时间轨道交会最优速度增量快速估计方法，适用于升交点赤经差较大而其它轨道根数较为接近的轨道交会问题。本方法可根据航天器与交会目标轨道根数差值和固定的转移时间，把复杂的多脉冲数值优化问题转化为简单等式/不等式约束的多元函数极值问题。本方法解决了无需精确求解时快速估计最优交会速度增量的问题，估计精度较高。
[0005] 为了实现上述技术目的，本发明所采用的技术方案是：
[0006] 一种J2摄动下长时间轨道交会最优速度增量快速估计方法，即已知航天器和目标在转移起始时刻t0的轨道根数的情形下，计算在tf时刻交会，转移时长为Δt的最优多脉冲速度增量Δvsum的过程。具体包括如下步骤：
[0007] (1)，计算需要通过脉冲消除的航天器和目标在交会时刻的轨道根数差；
[0008] (2)，设计消除ΔΩ0的最优控制模型，计算过脉冲消除ΔΩ0后交会所需的控制量以及最小速度增量；
[0009] (3)，将步骤(2)得到的解作为初值，计算高阶近似模型下交会所需的最小速度增量；
[0010] (4)，根据相位差对步骤(3)得到的最小速度增量进行修正，得到修正相位差后的交会所需的最小速度增量；
[0011] (5)，根据偏心率矢量差对步骤(4)得到最小速度增量进行修正，获得了估计出的最优交会速度增量。
[0012] 本发明步骤(1)中：已知航天器和目标在起始时刻t0的轨道根数。设tf为交会时刻，转移起始时刻到交会时刻之间的转移时长为Δt。其中轨道根数分别为轨道半长轴a、轨道偏心率e、轨道倾角i、升交点赤经Ω、近心点角ω、平近点角M。
[0013] 按照J2解析动力学模型分别计算交会时刻tf的航天器轨道根数和目标轨道根数。其中，J2解析动力学模型为：
[0014]
[0015] 其中： 为轨道根数对时间的导数，J2为地球引力二阶项，Re为地球半径，p＝a(1-e2)为轨道半通径， 为轨道角速度，μ为地球引力常数。
[0016] 则需要通过脉冲消除的航天器和目标在交会时刻的轨道根数差[Δa0,Δe0,Δi0,ΔΩ0,Δω0,ΔM0]为tf时刻的目标轨道根数减去航天器轨道根数。
[0017] 本发明步骤(2)中消除ΔΩ0的最优控制模型为：在转移起始时刻向航天器施加控制使航天器轨道根数产生变化量Δa,Δi,ΔΩ'，其中Δa为转移起始时刻航天器轨道半长轴的变化量，Δi为转移起始时刻航天器轨道倾角变化量，ΔΩ'为转移起始时刻航天器升交点赤经变化量；在交会时刻向航天器施加反向相同大小的控制量；Δa,Δi,ΔΩ'需要满足交会时刻航天器与目标相对升交点赤经为零的约束，在这个约束下求解最优的控制量中Δa,Δi,ΔΩ'分配，使速度增量取极小值。具体实现方法如下：
[0018] 对于航天器轨道半长轴a以及升交点赤经Ω，根据J2解析动力学模型即式(1)，在航天器与交会目标的轨道倾角接近，且轨道偏心率接近0的情形下，认为：
[0019]
[0020] 对航天器轨道倾角和轨道半长轴的控制量Δi和Δa会引起航天器升交点赤经的改变：
[0021]
[0022] 设脉冲施加方式为：在转移起始时刻tf-Δt施加单脉冲Δv0改变了航天器轨道半长轴和轨道倾角，进而改变航天器升交点赤经 然后在交会时刻tf施加单脉冲Δvf将航天器轨道半长轴和轨道倾角恢复为未施加控制量时的初始值，则tf时刻交会目标相对航天器的升交点赤经差为：
[0023]
[0024] 把消除ΔΩ′所需的速度增量一分为二合成到双脉冲里，则为了消除ΔΩ0，航天器需要的速度增量为：
[0025]
[0026] 式(5)中Δva，Δvi，ΔvΩ分别是航天器轨道半长轴变化量Δa，轨道倾角变化量Δi和升交点赤经变化量ΔΩ′对应的速度增量， 为航天器初始时刻轨道平均速度，且需要假设sini≈sin(i+Δi)。
[0027] 则最优速度增量的估算问题转为求Δvsum也即Δv0的极值问题。最优速度增量的估算问题可描述为：
[0028]
[0029] 将 的计算作一阶近似，则
[0030]
[0031] 其中 为航天器初始时刻升交点赤经漂移率，航天器初始时刻升交点赤经漂移率的计算方法为按照J2解析动力学模型即式(1)在航天器初始时刻轨道根数处线性展开只保留一阶项获得。
[0032] 考虑归一化变量x＝Δa/a，y＝Δi，z＝ΔΩ'，f和g等价于：
[0033]
[0034] 定义L＝f+λg，其中λ为约束乘子，则极值条件为：
[0035]
[0036] 式(9)是一个四阶线性方程，可以直接解得：
[0037]
[0038]
[0039] 则仅消除ΔΩ0所需的最小速度增量为
[0040]
[0041] 进一步地，本发明步骤(3)的实现方法如下：
[0042] 为了在tf时刻消除Δa0和Δi0，将式(5)中Δvf改写为：
[0043]
[0044] 同时考虑到Δa和Δi较大时，利用式(7)计算 的误差较大，采用式(3)计算并加入了交会机动过程中航天器轨道半长轴不能小于amin的约束(amin取为6600km)，得到更精确的求解模型为：
[0045]
[0046] 其中amin为约束的最小航天器轨道半长轴。
[0047] 同样设x＝Δa/ac，y＝Δi，z＝ΔΩ'，定义L＝f+λg1+χg2，其中λ和χ分别为约束乘子，相应的极值条件为：
[0048]
[0049] 其中 为式(2)中a替换为a+ax，i替换为i+y时计算的 即
[0050] 此时方程组(15)为非线性，将步骤(2)中公式(11)中求解到的x、y、z、λ作为方程组(15)中x、y、z、λ的初值，方程组(15)中χ初值设为0，进行多步迭代修正，得到新的x、y、z、λ和χ，进而得到新的
[0051]
[0052] 进一步地，本发明步骤(4)的实现方法如下：
[0053] 设Δλ0＝Δω0+ΔM0为无控制时航天器与目标在tf时刻的相位差，施加Δv0后，tf时刻相位差改变为：
[0054]
[0055] 将Δλ归一化到[-π,π]，为了满足交会约束，需要修正Δa和Δi产生一个相位的相对漂移率
[0056]
[0057] 考虑到sin2(i+Δi)对i的变化不敏感，认为Δa的变化起主导作用；则修正量为：
[0058]
[0059] 将其叠加到Δv0，即式(5)中固定xfix＝(Δa+Δamod)/a，方程组转化为：
[0060]
[0061] 使用步骤(3)中通过迭代修正后得到的新的y，z，λ作为式(20)中y，z，λ的初值，重新求解Δi和ΔΩ'，即可得到修正相位差后的交会所需的最小速度增量
[0062] 进一步地，本发明步骤(5)的实现方法如下：
[0063] 对航天器在初始时刻施加Δv0以后，航天器ω的漂移率有变化，对应交会时刻tf新的ω为：
[0064]
[0065] tf时刻航天器与目标的轨道偏心率差为
[0066]
[0067] 再根据Δv0、Δvf的法向分量对航天器轨道倾角i和升交点赤经Ω的改变量，可知这两个脉冲的控制相位为
[0068]
[0069] 其中在u0处施加Δv0对航天器轨道偏心率的改变量会随着ω的漂移一起旋转，认为tf时刻这个改变量对应的相位是：
[0070]
[0071] 如果将Δv0，Δvf分别拆分为间隔半个周期的双脉冲，则根据脉冲大小分配不同，其对航天器轨道偏心率的改变量Δev0和Δevf有：
[0072] -Δa/a＜|Δev0|＜Δa/a
[0073] -(Δa-Δa0)/a＜|Δevf|＜(Δa-Δa0)/a (25)
[0074] 设k0、kf∈[0,1]分别为Δv0，Δvf对航天器轨道偏心率的改变量系数，叠加后得Δv0，Δvf对航天器轨道偏心率的改变量Δe＝[Δex,Δey]为：
[0075]
[0076] 利用已有控制量来消除相对目标的偏心率差，则等同于求解极值问题：
[0077]
[0078] 解得min(f)以后，如果Δerem值大于零，就需要额外的轨道半长轴调整量：Δae＝aremΔe ，则按照式(28)在步骤(4)的Δvf中叠加即获得了估计出的最优交会速度增量[0079]
[0080] 相对于现有技术，本发明的有益技术效果是：
[0081] J2摄动主要带来升交点赤经Ω和近地点幅角ω的漂移，合理利用漂移可显著减小用于消除升交点赤经差和偏心率差所需的速度增量，这也是本发明的基本出发点。估计方法可应用于对多个目标连续交会顺序和时刻的全局优化，也可为多脉冲数值优化提供合理的初值，提高优化效率。
[0082] 本发明根据给定的航天器与交会目标轨道根数差值和转移时间，把复杂的多脉冲数值优化问题转化为简单等式/不等式约束的多元函数极值问题。解决了无需精确求解时快速估计最优交会速度增量的问题。
附图说明
[0083] 图1是本发明的流程图；
具体实施方式
[0084] 参照图1，为本发明的流程图，本发明提供一种J2摄动下长时间轨道交会最优速度增量快速估计方法。
[0085] 步骤1，计算航天器和目标在交会时刻的轨道根数差；
[0086] 已知任务航天器和目标在转移起始时刻t0的六个轨道根数分别为[ac(t0),ec(t0),ic(t0),Ωc(t0),ωc(t0),Mc(t0)]和[at(t0),et(t0),it(t0),Ωt(t0),ωt(t0),Mt(t0)]。
[0087] 其中ac(t0),ec(t0),ic(t0),Ωc(t0),ωc(t0),Mc(t0)分别代表航天器在转移起始时刻t0的六个轨道根数，分别为轨道半长轴、轨道偏心率、轨道倾角、升交点赤经、近心点角、平近点角：at(t0),et(t0),it(t0),Ωt(t0),ωt(t0),Mt(t0)分别代表目标在转移起始时刻t0的六个轨道根数，分别为轨道半长轴、轨道偏心率、轨道倾角、升交点赤经、近心点角、平近点角。
[0088] 设tf时刻为交会时刻，转移起始时刻到交会时刻之间转移时长为Δt；首先按照J2解析动力学模型分别获得tf时刻的航天器轨道根数
[0089] [ac(tf),ec(tf),ic(tf),Ωc(tf),ωc(tf),Mc(tf)]和目标轨道根数
[0090] [at(tf),et(tf),it(tf),Ωt(tf),ωt(tf),Mt(tf)]。
[0091] J2解析动力学模型如下：
[0092]
[0093] 其中， 为轨道根数对时间的导数，J2为地球引力二阶项，Re为地球半径，p＝a(1-e2)为轨道半通径， 为轨道角速度，μ为地球引力常数；
[0094] 则需要通过脉冲消除的航天器和目标在交会时刻的轨道根数差[Δa0,Δe0,Δi0,ΔΩ0,Δω0,ΔM0]为tf时刻的目标轨道根数减去航天器轨道根数。
[0095] 步骤2，计算通过脉冲仅消除ΔΩ0后交会所需的控制量；
[0096] 考虑仅消除ΔΩ0的最优控制模型为：在转移起始时刻向航天器施加控制使轨道根数产生变化量Δa,Δi,ΔΩ'，其中Δa为转移起始时刻航天器轨道半长轴的变化量，Δi为转移起始时刻航天器轨道倾角变化量，ΔΩ'为转移起始时刻航天器升交点赤经变化量。在交会时刻向航天器施加反向相同大小的控制量。Δa,Δi,ΔΩ'需要满足交会时刻航天器与目标相对升交点赤经为零的约束，在这个约束下求解最优的控制量中Δa,Δi,ΔΩ'分配，使速度增量取极小值。
[0097] 对于航天器轨道半长轴a以及航天器轨道升交点赤经Ω，根据式(1)，在航天器与交会目标的轨道倾角接近，且轨道偏心率接近0的情形下，可认为：
[0098]
[0099] 对航天器轨道倾角和轨道半长轴的控制量Δi和Δa会引起航天器升交点赤经的改变：
[0100]
[0101] 设脉冲施加方式为：在转移起始时刻tf-Δt施加单脉冲Δv0改变了航天器轨道半长轴和轨道倾角，进而改变航天器升交点赤经 然后在交会时刻tf施加单脉冲Δvf将航天器轨道半长轴和轨道倾角恢复为未施加控制量时的初始值，则tf时刻交会目标相对航天器的升交点赤经差为：
[0102]
[0103] 把消除ΔΩ′所需的速度增量一分为二合成到双脉冲里，则为了消除ΔΩ0，航天器需要的速度增量为：
[0104]
[0105] 式(5)中Δva，Δvi，ΔvΩ分别是航天器轨道半长轴变化量Δa，航天器轨道倾角变化量Δi和航天器升交点赤经变化量ΔΩ′对应的速度增量， 为航天器初始时刻轨道平均速度，且需要假设sini≈sin(i+Δi)。
[0106] 则最优速度增量的估算问题转为求Δvsum也即Δv0的极值问题。最优速度增量的估算问题可描述为：
[0107]
[0108] 将 的计算作一阶近似，则
[0109]
[0110] 其中 为航天器初始时刻升交点赤经漂移率。航天器初始时刻升交点赤经漂移率的计算方法为按照J2解析动力学模型即式(1)在航天器初始时刻轨道根数处线性展开只保留一阶项。
[0111] 考虑归一化变量x＝Δa/a，y＝Δi，z＝ΔΩ'，f和g等价于：
[0112]
[0113] 定义L＝f+λg(λ为约束乘子)，则极值条件为：
[0114]
[0115] 这是一个四阶线性方程，可以直接解得：
[0116]
[0117]
[0118] 则仅消除ΔΩ0所需的最小速度增量为
[0119]
[0120] 步骤3，计算同时消除[Δa0,Δi0,ΔΩ0]后交会所需的控制量。
[0121] 考虑同时消除Δa0,Δi0,ΔΩ0的最优控制模型为：在开始时刻施加控制量Δa,Δi,ΔΩ'，在交会时刻施加控制量为Δa0-Δa,Δi0-Δi,ΔΩ'。Δa,Δi,ΔΩ'需要满足交会时刻航天器与目标升交点赤经、半长轴和倾角均相等的约束，在这个约束下求解最优的Δa,Δi,ΔΩ'，使速度增量取极小值。
[0122] 由于还需要在tf时刻消除Δa0和Δi0，则式(5)中Δvf需要改写为：
[0123]
[0124] 同时考虑到Δa和Δi较大时，利用式(7)计算 的误差较大，采用式(3)计算并加入了交会机动过程中航天器轨道半长轴不能小于amin的约束(因为航天器地心距不能小于地球半径加上大气层高度，此处amin取为6600km)，得到更精确的求解模型为：
[0125]
[0126] 同样设x＝Δa/ac，y＝Δi，z＝ΔΩ'，定义L＝f+λg1+χg2(λ和χ分别为约束乘子)，相应的极值条件为：
[0127]
[0128] 其中 为式(2)中a替换为a+ax，i替换为i+y时计算的 即
[0129] 此时方程组(15)为非线性，通过步骤2中公式(11)中求解到的x、y、z、λ作为方程组(15)中x、y、z、λ的初值，方程组(15)中χ初值设为0，进行多步迭代修正，得到新的x、y、z、λ和χ，进而得到新的
[0130]
[0131] 步骤4，根据相位差对步骤3中得到的控制量进行修正。
[0132] 设Δλ0＝Δω0+ΔM0为无控制时航天器与目标在tf时刻的相位差，施加Δv0后，tf时刻相位差改变为：
[0133]
[0134] 将Δλ归一化到[-π,π]，为了满足交会约束，需要修正Δa和Δi产生一个相位的相对漂移率
[0135]
[0136] 考虑到sin2(i+Δi)对i的变化不敏感，可认为Δa的变化起主导作用。则修正量为：
[0137]
[0138] 将其叠加到Δv0，即式(5)中固定xfix＝(Δa+Δamod)/a，方程组转化为：
[0139]
[0140] 使用步骤3中通过迭代修正后得到的新的y，z，λ作为式(20)中y，z，λ的初值，重新求解Δi和ΔΩ'，即可得到修正相位差后的交会所需的最小速度增量
[0141] 步骤5，根据偏心率矢量差对步骤4中得到的控制量进行修正。
[0142] 对航天器在初始时刻施加Δv0以后，航天器ω的漂移率有变化，对应交会时刻tf新的ω为：
[0143]
[0144] tf时刻航天器与目标的轨道偏心率差为
[0145]
[0146] 再根据Δv0、Δvf的法向分量对航天器轨道倾角i和升交点赤经Ω的改变量，可知这两个脉冲的控制相位为：
[0147]
[0148] 其中在u0处施加Δv0对航天器轨道偏心率的改变量会随着ω的漂移一起旋转，所以可以认为tf时刻这个改变量对应的相位是：
[0149]
[0150] 如果将Δv0，Δvf分别拆分为间隔半个周期的双脉冲，则根据脉冲大小分配不同，其对航天器轨道偏心率的改变量Δev0和Δevf有：
[0151] -Δa/a＜|Δev0|＜Δa/a
[0152] -(Δa-Δa0)/a＜|Δevf|＜(Δa-Δa0)/a (25)
[0153] 设k0、kf∈[0,1]分别为Δv0，Δvf对航天器轨道偏心率的改变量系数，叠加后得Δv0，Δvf对航天器轨道偏心率的改变量Δe＝[Δex,Δey]为
[0154]
[0155] 要尽量利用已有控制量来消除相对目标的偏心率差，则等同于求解极值问题：
[0156]
[0157] 解得min(f)以后，如果Δerem值大于零，就需要额外的轨道半长轴调整量：Δae＝aΔerem，则按照式(28)在步骤(4)的Δvf中叠加即获得了估计出的最优交会速度增量[0158]
[0159] 下面结合一具体实例，采用本发明方法获取最优速度增量。
[0160] 本实施例中：航天器和目标初始时刻t0＝0轨道根数见表1：
[0161] 表1航天器和目标初始时刻轨道根数
[0162] A(m) e i(rad) Ω(rad) w(rad) M(rad)
航天器 6996104.270 0.0016195 1.72108 3.79428 0.27199 3.09317
目标 7086789.952 0.005314 1.73557 1.61260 4.46103 0.49992
[0163] 要计算在第2365天开始轨道转移，第2390天交会所需的最优速度增量。即Δt＝25天，tf为第2390天。
[0164] 第一步，计算航天器和目标在交会时刻的轨道根数差。
[0165] 首先将轨道根数外推至tf时刻，得轨道根数差[Δa0,Δi0,ΔΩ0]为[90685.682m 0.01449 -0.0501]。
[0166] 第二步，根据一阶近似模型计算交会所需的控制量。
[0167] 根据ΔΩ0和式(11)求解得到[x y z]＝[0.0153054 -0.007219 -0.002374]。则一阶近似模型下消除升交点赤经差所需的初始时刻最优控制量[Δa,Δi,ΔΩ']分别为[107078.214m -0.007219 -0.002374]。
[0168] 第三步，将第二步的解作为初值，计算高阶近似模型下交会所需的控制量。
[0169] 将其作为初值，代入式(15)，其中amin设为地球半径加上200km，解得[x y z]＝[0.02123 -0.001943 -0.01105]。代入式(14)，得消除[Δa0,Δi0,ΔΩ0]所需的最优控制量Δv0和Δvf分别为90.4m/s和133.1m/s。
[0170] 第四步，根据相位差对第三步得到的控制量进行修正。
[0171] 根据航天器轨道根数和Δv0，可重新计算在tf时刻目标相对航天器相位差为1.8419rad，对应Δt所需的半长轴修正量为-3686.3m。将其代入xfix＝(Δa+Δamod)/a＝
0.0207，求解式(20)，得[y z]＝[-0.002154 -0.01122]。重新代入式(14)，即修正相位差后最优控制量Δv0和Δvf分别为89.3m/s和134.3m/s。
[0172] 第五步，根据偏心率矢量差对第四步得到的控制量进行修正。
[0173] 根据航天器轨道根数和Δv0，计算得在tf时刻目标相对航天器偏心率差为：由式(23)得u0和uf分别为-2.6762和-1.8925，求解式(27)得
Δerem＝0，即Δae＝aΔerem为0。
[0174] 因此估计出的最优交会速度增量为Δvsum＝Δv0+Δvf＝223.6m/s。使用J2摄动动力学模型按照四脉冲数值优化计算得到得速度增量为226.8m/s。
[0175] 经实例论证，采用本发明方法估计速度增量的误差较小。本方法解决了无需精确求解时快速估计最优交会速度增量的问题，估计精度较高。
[0176] 以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。
