具体实施方式

下面将参照附图更详细地描述本公开的示例性实施例。虽然附图中显示了本公开的示例性实施例，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。

本实施例提供了一种化电混合推进变轨优化方法，先化推后电推，化推段变轨按照脉冲变轨将变轨点位置、变轨推力方向、变轨速度增量作为优化变量；电推段采用近似等效方法，仅对2个控制参数进行优化：轨控仰角Ψ₁、Ψ₂，且可以通过半解析方法求解。最终将化电混合推进变轨按照参数优化问题进行求解，优化参数不超过10个，计算过程简单，计算结果可靠，变轨指标优于常规设计值。

本发明将变轨过程分为化推变轨、电推变轨2个阶段。首先化推段按照双脉冲变轨等效，从初始轨道变轨至中间轨道，变轨参数包括：变轨点真近点角θ₁、θ₂，变轨推力仰角变轨速度增量ΔV₁、ΔV₂。电推段按照近似等效方法，以化推段变轨结束的中间轨道作为起始轨道，变轨参数包括：轨控仰角Ψ₁、Ψ₂，可以通过半解析方法求解最优参数Ψ₁ ^*、Ψ₂ ^*。化电混合变轨可以根据6个变轨参数θ₁、θ₂、ΔV₁、ΔV₂计算得到变轨时间或推进剂消耗指标，采用优化算法得到最优参数θ₁ ^*、θ₂ ^*、ΔV₁ ^*、ΔV₂ ^*，再结合电推段变轨参数Ψ₁ ^*、Ψ₂ ^*得到化电混合推进最优变轨方法。

该方法具体包括以下步骤：

(1)步骤一、确定变轨过程初末状态参数。变轨过程分为化推变轨、电推变轨2个阶段，先进行化推变轨、后进行电推变轨。采用经典轨道根数a、e、i、Ω、ω、θ(半长轴、偏心率、倾角、升交点赤经、近地点幅角、真近点角)，其中真近点角仅表征卫星相位，描述卫星初末轨道时不予考虑，初始轨道参数为(a₀,e₀,i₀,Ω₀,ω₀)、目标轨道参数为(a_f,e_f,i_f,Ω_f,ω_f,θ_f)，中间轨道即化推电推变轨分割点参数记为(a_m,e_m,i_m,Ω_m,ω_m)为待确定参数。

(2)步骤二、确定化推变轨方程。化学推进大推力作用下，变轨时间一般较短，点火变轨的弧段损失可忽略，可等效为脉冲变轨。化推段变轨采用双脉冲变轨，第1个脉冲在近地点附近，速度增量为ΔV₁，点火点的真近点角为θ₁，脉冲方向沿轨道切向，并向轨道法向抬起仰角γ₁；第2个脉冲在远地点附近，速度增量为ΔV₂，点火点的真近点角为θ₂，脉冲方向沿轨道切向，并向轨道法向抬起仰角γ₂。

第1次脉冲点火前卫星轨道参数为(a₀,e₀,i₀,Ω₀,ω₀,θ₁)，转换成位置速度矢量：r₁₀、V₁₀，对应的轨道法向矢量n₁，点火后卫星位置速度矢量变为r_1f、V_1f：

点火后位置速度矢量r_1f、V_1f再转换为点火后轨道参数记为(a_k,e_k,i_k,Ω_k,ω_k,θ_k)。

然后，第2次脉冲点火前卫星轨道参数为(a_k,e_k,i_k,Ω_k,ω_k,θ₂)，转换成位置速度矢量：r₂₀、V₂₀，对应的轨道法向矢量n₂，点火后卫星位置速度矢量变为r_2f、V_2f：

点火后位置速度矢量r_2f、V_2f再转换为点火后轨道参数即为中间轨道参数(a_m,e_m,i_m,Ω_m,ω_m,θ_m)。

综上，化推变轨方程轨道参数变化与变轨参数θ₁、θ₂、ΔV₁、ΔV₂的关系，可记为：

化推变轨方程(a_m,e_m,i_m,Ω_m,ω_m)＝f(a₀,e₀,i₀,Ω₀,ω₀,θ₁,θ₂,γ₁,γ₂,ΔV₁,ΔV₂)

其中，真近点角对电推变轨阶段计算无影响，遂略去。

(3)步骤三、确定电推变轨方程。电推变轨过程分为3个阶段：

a)第1阶段，推力沿轨道切向，将近地点高度抬高至1000km，以减小大气阻力摄动影响。

b)第2阶段，推力沿轨道切向，并向轨道法向抬起轨控仰角Ψ₁，在轨道幅角为90°和270°前后，Ψ₁改变正负，保证推力方向在升交点前后偏向轨道负法向、在降交点前后偏向轨道正法向，将半长轴增大至42164km。

c)第3阶段，推力方向沿过近地点速度的反方向并向轨道法向抬起轨控仰角Ψ₂，在轨道幅角为90°和270°前后，Ψ₁改变正负，保证推力方向在升交点前后偏向轨道负法向、在降交点前后偏向轨道正法向，将偏心率降至0，同时倾角也降至约0°。

分阶段给出电推变轨动力学方程。第1阶段变轨动力学方程：

其中n为轨道平均速率，F_p为变轨总推力，m为卫星质量，J₂为地球非球形摄动二阶项常数，R_e为地球半径，k_a ^I、k_e ^I为第1阶段变轨过程半长轴a、偏心率e的轨道平均因子，只与当前轨道的偏心率e有关：

将轨道根数变化转换为以偏心率为自变量：

采用第1阶段变轨动力学方程进行轨道积分计算，得到第1阶段变轨速度增量ΔV^I：

其中e₁为第1阶段变轨过程结束时刻的偏心率。

第2阶段变轨轨迹、变轨指标与待确定变轨参数——轨控仰角Ψ₁有关，因此不能直接采用变轨动力学方程的积分方式得到变轨轨迹、变轨时间。需要先计算变轨第2阶段重要积分系数A。第2阶段变轨动力学方程为：

其中k_a ^II、k_e ^II、k_i ^II、k_ω ^II为第2阶段变轨过程半长轴a、偏心率e、倾角i、近地点幅角ω的轨道平均因子，k_a ^II、k_e ^II只与当前轨道的偏心率e有关，k_i ^II、k_ω ^II与当前轨道的偏心率e、近地点幅角ω有关：

将轨道根数变化转换为以偏心率为自变量：

采用上式对第2阶段轨道参数的变化进行计算，并计算第2阶段重要积分系数A，以及第2阶段变轨速度增量在轨道平面内的分量

其中e₂为第2阶段变轨过程结束时刻的偏心率。

同理，第3阶段变轨轨迹、变轨指标与待确定的变轨参数：轨控仰角Ψ₂有关，因此不能直接采用变轨动力学方程的积分方式得到变轨轨迹、变轨时间。需要先计算变轨第3阶段重要积分系数B。第3阶段变轨动力学方程：

其中k_e ^III、k_i ^III、k_ω ^III为第3阶段变轨过程偏心率e、倾角i、近地点幅角ω的轨道平均因子，k_e ^III只与当前轨道的偏心率e有关，k_i ^III、k_ω ^III与当前轨道的偏心率e、近地点幅角ω有关：

将轨道根数变化转换为以偏心率为自变量：

采用上式对第3阶段轨道参数的变化进行积分计算，并计算第3阶段重要积分系数B，以及第3阶段变轨速度增量在轨道平面内的分量

最优变轨参数Ψ₁ ^*、Ψ₂ ^*为如下方程组的解：

求解上式得到最优变轨参数Ψ₁ ^*、Ψ₂ ^*。

由于近地点幅角、倾角变化的耦合作用较为突出，而近地点幅角受到地球非球形引力J₂项摄动影响较为明显，在计算过程中需要迭代2～3次才可得到精确的计算结果。

综上，电推变轨方程描述轨道参数变化及变轨速度增量ΔV_E与最优变轨参数Ψ₁ ^*、Ψ₂ ^*的关系，可记为：

(4)步骤四、确定化电混合变轨优化模型并求解最优变轨参数。

根据上述变轨过程描述，化推变轨后进入中间轨道，后续电推变轨段最优变轨参数Ψ₁ ^*、Ψ₂ ^*及对应的速度增量ΔV_E可以采用半解析方法求解，即整个化电混合变轨的优化只需要对化推变轨参数θ₁、θ₂、ΔV₁、ΔV₂进行优化，优化模型为：

其中：J为优化指标，化推变轨阶段速度增量与电推变轨阶段速度增量的加权和；λ为化电变轨优化的权重系数，取值0～1，λ＝0时，即要求在给定化推推进剂消耗的条件下使得变轨时间最短，λ＝1时，即要求在给定变轨时间条件下使得推进剂消耗最省；ΔV_U ^upp、ΔV_E ^upp分别为化推变轨、电推变轨可用的速度增量上限约束。Ψ₁ ^*、Ψ₂ ^*与中间轨道参数有关，即需满足电推变轨方程：

中间轨道参数与化推段变轨参数有关，即需满足化推变轨方程：

(a_m,e_m,i_m,Ω_m,ω_m)＝f(a₀,e₀,i₀,Ω₀,ω₀,θ₁,θ₂,γ₁,γ₂,ΔV₁,ΔV₂) (18)

采用优化算法对式(16)进行求解，即可得到化推段最优变轨参数θ₁ ^*、θ₂ ^*、ΔV₁ ^*、ΔV₂ ^*。

(5)步骤五、根据最优变轨参数计算化推、电推变轨轨迹及相关参数。

根据步骤四得到的最优变轨参数，再次运用化推变轨方程、电推变轨方程，计算得到中间轨道参数(a_m,e_m,i_m,Ω_m,ω_m,θ_m)、电推变轨段速度增量ΔV_E，再利用火箭方程，可计算得到化推推进剂消耗Δm_U,电推推进剂消耗Δm_E以及电推变轨时间ΔT_E：

其中，m₀为卫星初始重量，I_sp ^U、I_sp ^E分别为化推、电推比冲，单位为m/s。

由上式可见，只针对电推变轨段而言，变轨时间ΔT_E与速度增量ΔV_E唯一相关，因此，步骤四中优化模型的优化指标、约束函数可以根据实际任务需求，将ΔV_E替换为ΔT_E。

GEO卫星发射重量5400kg，初始轨道Hp＝200km，Ha＝35786km，i＝28.5°，w＝179°。化学推进等效为脉冲变轨，比冲312s；电推进推力500mN，比冲3000s。

设计两组化电混合推进变轨任务算例分别为：

1)算例一：以化推变轨推进剂消耗最省为目标，限定电推变轨时间上限180天。

2)算例二：以电推进变轨时间最短为目标，限定化推推进剂消耗上限1500kg。

(1)步骤一、确定变轨过程初末状态参数。初始轨道参数(a₀,e₀,i₀,Ω₀,ω₀)＝(24371km,0.7301,28.5°,0°,179°)，目标轨道参数(a_f,e_f,i_f,Ω_f,ω_f)＝(42164km,0,0°,0°,0°)。中间轨道参数(a_m,e_m,i_m,Ω_m,ω_m)待定。

(2)步骤二、确定化推变轨方程。

(3)步骤三、确定电推变轨方程。

(4)步骤四、确定化电混合变轨优化模型并求解最优变轨参数。对于算例一，优化目标为化推变轨推进剂消耗最省，约束条件为电推变轨时间上限180天，令λ＝1，并将ΔV_E替换为ΔT_E，优化模型如下：

对于算例二，优化目标为电推进变轨时间最短，约束条件为化推推进剂消耗上限1500kg，转换成速度增量为995.6m/s，令λ＝0，并将ΔV_E替换为ΔT_E，优化模型如下：

(4)步骤五、采用优化算法求解，得到最优变轨参数及变轨过程相关参数如下表。

半长轴、偏心率、倾角的时间历程如图1(a)、图1(b)、图1(c)、图3(a)、图3(b)、图3(c)所示，图2、图4为变轨过程三维飞行轨迹。

本发明先化推后电推，化推段变轨按照脉冲变轨将变轨点位置、变轨推力方向、变轨速度增量作为优化变量；电推段采用近似等效方法，仅对2个控制参数进行优化：轨控仰角Ψ₁、Ψ₂，且可以通过半解析方法求解。最终将化电混合推进变轨按照参数优化问题进行求解，优化参数仅6个，计算过程简单，计算结果可靠，变轨指标优于常规设计值；本发明通过权重调节可以对化电混合推进变轨任务的分配进行灵活调整，适应多种任务需求，应用场景广泛。

本发明虽然已以较佳实施例公开如上，但其并不是用来限定本发明，任何本领域技术人员在不脱离本发明的精神和范围内，都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改，因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰，均属于本发明技术方案的保护范围。