CN105373133A - 一种同步轨道电推进位置保持与角动量卸载联合控制方法 - Google Patents

Abstract

本发明涉及一种同步轨道电推进位置保持与角动量卸载联合控制方法，包括：确定位置保持周期T_s、控制周期T_c和测轨时间T_m；控制周期T_c为包含所有推力器各完成一次点火过程的总时间；测轨时间T_m由测定轨所需时间向上取整为整数个轨道周期得到；一个位置保持周期T_s应包含一个测轨过程T_m和若干个控制周期T_c：每个位置保持周期T_s开始之前，计算每个控制周期轨道要素控制量：倾角控制量矢量(Δi_x,Δi_y)、偏心率控制量矢量(Δe_x,Δe_y)和平经度漂移率控制量ΔD；确定每个控制周期的电推力器点火参数，点火参数包括点火时长及点火位置；每个控制周期T_c开始之前，根据电推力器点火参数及角动量卸载量确定电推力器转动角度。

Description

一种同步轨道电推进位置保持与角动量卸载联合控制方法

技术领域

本发明涉及一种同步轨道卫星电推进位置保持与动量轮卸载联合控制方法，属于卫星轨道姿态动力学与控制技术领域。

背景技术

地球同步轨道(GEO，GeosynchronousOrbit)是指卫星运行时相对地球位置基本不变，与地球自转运动同步的卫星轨道。实际上由于各种摄动因素的作用，同步轨道卫星并不能始终保证相对地球位置不变，需要不断地施加控制力使卫星相对地球位置保持在要求的范围内。同步轨道卫星的位置保持操作实际上是对卫星倾角、偏心率和平经度的控制。

目前化学推进卫星的位置保持是通过推力方向解耦的化学推力器根据需要分别产生沿轨道径向、切向、法向的速度增量对同步轨道卫星的定点位置进行控制的。化学推进卫星的角动量卸载则是通过推力器组合使其对星体的推力相互抵消，只产生相应的卸载力矩实现的。其位置保持与角动量卸载相互解耦。

目前采用电推进技术的同步轨道卫星，其电推力器布局一般使其推力在轨道径向、切向、法向的分量均不为零，即在南北、东西方向的位置保持控制上存在耦合；且其角动量卸载也需要电推力器执行，电推力器在产生绕星体质心的力矩的同时必定会产生对星体的推力从而引起轨道参数变化，即位置保持控制与角动量卸载存在耦合。解耦控制会使得推进剂利用效率低，造成推进剂的浪费。

北京理工大学CN104015938A号专利提供了一种电推进静止轨道卫星的位置保持方法，该方法通过将一个位置保持周期划分为若干个小控制周期，每个位置保持周期测轨数据更新一次，根据测轨数据计算每个小控制周期的电推力器点火参数，实现东西、南北方向位置保持的联合控制，按照给定的位置保持周期重复执行，达到了较高的位置保持精度和较高的推进剂利用效率。但是该方法主要存在以下问题：

(1)未考虑位置保持控制参数(位置保持周期、控制周期等)与位置保持精度之间的关系，在控制参数的选择上存在盲目性。

(2)在每个位置保持周期内，不同的控制周期的电推力器点火参数需要递推计算求得，使得各控制周期的点火参数并不相同，解算复杂。

(3)未考虑角动量卸载需求，未考虑电推力器执行角动量卸载对位置保持的影响。

发明内容

本发明的技术解决问题是：克服以上缺点，提供一种简单易行的同步轨道电推进位置保持与角动量卸载联合控制方法，解决同步轨道电推进位置保持与角动量卸载问题，并实现较高的推进剂利用效率。

本发明的技术解决方案是：一种同步轨道电推进位置保持与角动量卸载联合控制方法，以星上测轨数据的更新周期为一个位置保持周期T_s，以全部电推力器各点火一次为一个控制周期T_c，一个位置保持周期包含若干个控制周期。首先在一个置保持周期开始之前根据测轨数据确定一个位置保持周期所需的轨道要素控制量，将其平均分配到每个控制周期，确定每个控制周期内每台电推力器的点火时长以及点火位置；其次在每个控制周期开始之前，根据动量轮系的角动量卸载需求，调整该控制周期内每台电推力器的指向，以产生相应的卸载力矩。按照位置保持周期重复执行该策略，即可实现卫星定点位置和轮系角动量的控制。

实现步骤如下：

步骤一、确定位置保持周期T_s、控制周期T_c和测轨时间T_m

控制周期T_c为包含所有推力器各完成一次点火过程的总时间，可取1天或2天；测轨时间T_m由测定轨所需时间向上取整为整数个轨道周期得到；一个位置保持周期T_s应包含一个测轨过程T_m和若干个控制周期T_c：

T_s＝T_m+mT_c

其中m为正整数；

步骤二、每个位置保持周期T_s开始之前，根据测轨得到的轨道数据计算每个控制周期轨道要素控制量：倾角控制量矢量(Δi_x,Δi_y)、偏心率控制量矢量(Δe_x,Δe_y)和平经度漂移率控制量ΔD；

步骤三、根据各轨道要素所需控制量确定每个控制周期的电推力器点火参数，点火参数包括点火时长及点火位置；

步骤四、每个控制周期T_c开始之前，根据电推力器点火参数及角动量卸载量确定电推力器转动角度。

重复步骤二到步骤四，在电推力器控制下卫星将保持在定点位置附近，且每个控制周期内角动量能够持续进行卸载，不存在角动量长期积累。

所述步骤一中m由东西方向位置保持精度要求确定，东西方向位置保持精度α_EW由偏心率控制精度e_d和平经度控制精度λ_d决定，如下：

α_EW＝2e_d(T_m)+λ_d(T_m,m)

m由上式求根并向下取整得到。

所述偏心率控制精度e_d按下式计算：

e_d＝Δe_meas+Δe_ctrl+Δe_rand+Δe_pert

其中Δe_meas、Δe_ctrl、Δe_rand、Δe_pert分别为偏心率测量误差、偏心率控制误差、不确定误差、偏心率漂移量；不确定误差Δe_rand包含了控制律中未予考虑的其他摄动项的影响，主要是日周期项，该项约为0.9×10^-4；偏心率漂移量Δe_pert表示无控弧段偏心率自由摄动所能达到的最大值，该项约为：

${\mathrm{\Δe}}_{pert}=max(T_m,\frac{1}{2}day)n_s{R}_e$

其中n_s为太阳视运动角速率，R_e＝0.011C_rS/m为偏心率摄动圆半径，C_r为光压系数，S/m为卫星面积质量比。

所述平经度控制精度λ_d按下式计算：

${\mathrm{\λ}}_d={\mathrm{\λ}}_{meas}+{\mathrm{\λ}}_{ctrl}+\frac{3{\mathrm{\πK}}_T}{4K_N}{\mathrm{m\Δi}}_C+\frac{{\stackrel{\·\·}{\mathrm{\λ}}}_s{T}_m{T}_s}{16}$

其中λ_meas、λ_ctrl分别为平经度测量误差、平经度控制误差；Δi_C为一个控制周期内的倾角控制量，当控制周期T_c取为1天，Δi_C为0.0026°，当控制周期T_c取为2天，Δi_C为0.0052°；ω_e为地球自转角速率，为定点位置处平经度漂移加速度；K_T、K_N分别为电推力器推力在轨道切向、法向的投影系数，取正值。

所述步骤二中的倾角控制量(Δi_x,Δi_y)按照如下控制律计算：

其中为太阳平黄经，t为距春分时间；A_x＝0.0235°，A_y＝0.0215°分别为倾角矢量半年周期摄动项的幅值；Ω_ms为白道升交点黄经，其变化可记为Ω_ms＝12.11°-0.052954T，其中T为相对儒略日，相对于1950年1月1日。

所述步骤二中偏心率控制量(Δe_x,Δe_y)按照如下控制律计算：

其中i_s＝23.45°为黄赤交角，为太阳平黄经，t为距春分点时间；R_e＝0.011C_rS/m为偏心率摄动圆半径，C_r为光压系数，S/m为卫星面积质量比；e_x0，e_y0为偏心率矢量测量值。

所述步骤二中和平经度漂移率控制量ΔD，一个位保周期内m个控制周期分成两段，前m/2个控制周期产生相同的平经度漂移率控制量ΔD₁，后m/2个控制周期产生相同的平经度漂移率控制量ΔD₂，按照如下步骤计算：

(1)计算平经度控制目标。平经度控制目标为在当前位保周期结束时平经度和平经度漂移率控制到指定值(λ^*,D^*)，(λ^*,D^*)可由下式计算：

${\mathrm{\λ}}^{*}=-\frac{5}{4}{\stackrel{\·\·}{\mathrm{\λ}}}_s\·{\left(1day\right)}^2$

$D^{*}=-{\stackrel{\·\·}{\mathrm{\λ}}}_s\·1day$

其中为定点位置处平经度漂移加速度。

(2)求解如下无约束优化问题：

$\min J=\min \left|\right|\left[\begin{array}{cc}\left(\frac{3m^2}{8}+\frac{m}{4}\right)T_c& \left(\frac{m^2}{8}+\frac{m}{4}\right)T_c\\ \frac{m}{2}{T}_c& \frac{m}{2}{T}_c\end{array}\right]X+\left[\begin{array}{c}{\mathrm{\λ}}_0+D_0{\mathrm{mT}}_c+\frac{1}{2}{\stackrel{\·\·}{\mathrm{\λ}}}_s{\left({\mathrm{mT}}_c\right)}^2-{\mathrm{\λ}}^{*}\\ D_0{T}_c+{\stackrel{\·\·}{\mathrm{\λ}}}_s{\mathrm{mT}}_c^2-D^{*}{T}_c\end{array}\right]|{|}_2$

X＝(ΔD₁,ΔD₂)^T

其中λ₀、D₀分别为平经度、平经度漂移率的测量值。求得一组最优解X^*＝(ΔD₁,ΔD₂)^T作为平经度漂移率的控制量。

所述步骤三中电推力器点火参数计算过程如下：

(1)首先根据倾角控制量(Δi_x,Δi_y)计算各电推力器点火弧段中点赤经：

l₁＝atan2(-Δi_y,-Δi_x)

l₂＝atan2(-Δi_y,-Δi_x)

l₃＝atan2(-Δi_y,-Δi_x)+π

l₄＝atan2(-Δi_y,-Δi_x)+π

(2)计算各电推力器所需产生的速度增量：

${\mathrm{\ΔV}}_1=\frac{1}{4}\left(\frac{V_s\mathrm{\Δ}i}{K_N}+\frac{V_s\left({\mathrm{\Δe}}_x{\mathrm{cosl}}_1+{\mathrm{\Δe}}_y{\mathrm{sinl}}_1\right)}{2K_T}+\frac{V_s\left(-{\mathrm{\Δe}}_x{\mathrm{sinl}}_1+{\mathrm{\Δe}}_y{\mathrm{cosl}}_1\right)}{K_R}-\frac{R_s\mathrm{\Δ}D}{3K_T}\right)$

${\mathrm{\ΔV}}_2=\frac{1}{4}\left(\frac{V_s\mathrm{\Δ}i}{K_N}-\frac{V_s\left({\mathrm{\Δe}}_x{\mathrm{cosl}}_1+{\mathrm{\Δe}}_y{\mathrm{sinl}}_1\right)}{2K_T}+\frac{V_s\left(-{\mathrm{\Δe}}_x{\mathrm{sinl}}_1+{\mathrm{\Δe}}_y{\mathrm{cosl}}_1\right)}{K_R}+\frac{R_s\mathrm{\Δ}D}{3K_T}\right)$

${\mathrm{\ΔV}}_3=\frac{1}{4}\left(\frac{V_s\mathrm{\Δ}i}{K_N}-\frac{V_s\left({\mathrm{\Δe}}_x{\mathrm{cosl}}_1+{\mathrm{\Δe}}_y{\mathrm{sinl}}_1\right)}{2K_T}-\frac{V_s\left(-{\mathrm{\Δe}}_x{\mathrm{sinl}}_1+{\mathrm{\Δe}}_y{\mathrm{cosl}}_1\right)}{K_R}-\frac{R_s\mathrm{\Δ}D}{3K_T}\right)$

${\mathrm{\ΔV}}_4=\frac{1}{4}\left(\frac{V_s\mathrm{\Δ}i}{K_N}+\frac{V_s\left({\mathrm{\Δe}}_x{\mathrm{cosl}}_1+{\mathrm{\Δe}}_y{\mathrm{sinl}}_1\right)}{2K_T}-\frac{V_s\left(-{\mathrm{\Δe}}_x{\mathrm{sinl}}_1+{\mathrm{\Δe}}_y{\mathrm{cosl}}_1\right)}{K_R}+\frac{R_s\mathrm{\Δ}D}{3K_T}\right)$

其中K_T、K_R、K_N分别为卫星对地姿态下电推力器推力沿轨道切向、径向、法向的投影系数，均取正值。

(3)根据各电推力器所需产生的速度增量计算实际的电推力器点火时长：

${\mathrm{\Δt}}_i=\frac{2}{{\mathrm{\ω}}_e}arcsin\left(\frac{m_c{\mathrm{\ω}}_e{\mathrm{\ΔV}}_i}{2F_p}\right)$

其中m_c为航天器重量，F_p为电推力器推力。

(4)根据电推力器点火时长与电推力器点火弧段中点赤经计算电推力器起始点火时刻：

$t_{start}=t_0+\frac{l_i}{{\mathrm{\ω}}_e}-\frac{1}{2}{\mathrm{\Δt}}_i$

其中t₀为过春分点时刻，ω_e为地球自转角速率。

所述步骤四中根据电推力器点火参数及角动量卸载量确定电推力器转动角度，计算过程如下：

(1)求解如下二次规划问题：

minJ＝X^TX

s.t.ΔH_C＝AX

X＝(Δθ₁,Δθ₂,Δθ₃,Δθ₄,Δα₁,Δα₂,Δα₃,Δα₄)^T

其中A∈R^3×8，A的第k列为：

$A_k=\left\{\begin{array}{c}{m}_{c}d\underset{i=1}{\overset{4}{\Σ}}{\mathrm{\ΔV}}_i{C}_x\left(\frac{\mathrm{\π}}{2}\right)C_y\left(l_i+\frac{\mathrm{\π}}{2}\right)\left[\begin{array}{c}-{\mathrm{cos\α}}_i\\ 0\\ {\mathrm{sin\α}}_i\end{array}\right],k=1,2,3,4\\ m_{c}d\underset{i=1}{\overset{4}{\Σ}}{\mathrm{\ΔV}}_i{C}_x\left(\frac{\mathrm{\π}}{2}\right)C_y\left(l_i+\frac{\mathrm{\π}}{2}\right)\left[\begin{array}{c}{\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i{\mathrm{sin\α}}_i\\ {\sin }^2{\mathrm{\θ}}_i\\ {\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i{\mathrm{cos\α}}_i\end{array}\right],k=5,6,7,8\end{array}\right.$

其中θ为推力方向与Y轴的夹角；α为推力在XOZ面的投影与Z轴的夹角，α_i、θ_i分别为各电推力器初始安装角度；d为电推力器作用点到质心的距离。C_x为绕坐标系X轴旋转的方向余弦矩阵，C_y为绕坐标系Y轴旋转的方向余弦矩阵，α_i、θ_i分别为各电推力器初始安装角度。解得一组最优解X⁰＝(Δθ₁,Δθ₂,Δθ₃,Δθ₄,Δα₁,Δα₂,Δα₃,Δα₄)^T为电推力器转角的初值。

(2)根据电推力器转角初值X⁰计算角动量卸载截断误差ΔH_2d：

${\mathrm{\ΔH}}_{2d}=m_{c}d\Σ{\mathrm{\ΔV}}_i{C}_x\left(\frac{\mathrm{\π}}{2}\right)C_y\left(l_i+\frac{\mathrm{\π}}{2}\right)\left(-\frac{1}{2}\left[\begin{array}{c}-{\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i{\mathrm{cos\α}}_i\\ 0\\ {\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i{\mathrm{sin\α}}_i\end{array}\right]{\mathrm{\Δ\α}}_i^2+\left[\begin{array}{c}{\cos }^2{\mathrm{\θ}}_i{\mathrm{sin\α}}_i\\ -{\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i\\ {\cos }^2{\mathrm{\θ}}_i{\mathrm{cos\α}}_i\end{array}\right]{\mathrm{\Δ\θ}}_i{\mathrm{\Δ\α}}_i\right)$

(3)考虑角动量卸载截断误差ΔH_2d重新求解如下二次规划问题：

minJ＝X^TX

s.t.ΔH_C+ΔH_2d＝AX

X＝(Δα₁,Δα₂,Δα₃,Δα₄,Δθ₁,Δθ₂,Δθ₃,Δθ₄)^T

得到一组最优解X^*＝(Δθ₁,Δθ₂,Δθ₃,Δθ₄,Δα₁,Δα₂,Δα₃,Δα₄)^T作为最终的电推力器转角值。

本发明与现有技术相比的优点在于：

(1)提供一种同步轨道电推进位置保持与角动量卸载联合控制方法，能够同时实现位置保持高精度控制和角动量管理，提高了推进剂利用效率，节省推进剂消耗量；

(2)本发明所提供的位置保持方法，能够一次性确定一个位保周期内的全部点火参数，计算简单；本发明所提供的角动量卸载方法，在每个控制周期之前，依据星上角动量卸载需求，灵活调整推力器转角，既能实现角动量的高精度控制，同时不影响位保策略的执行。

(3)本发明同时提供了位置保持周期、控制周期、测轨时间的确定方法，避免了参数选择的盲目性。

(4)本发明所提供的同步轨道电推进位置保持与角动量卸载联合控制方法，实施简单，自主性强，在具备自主测轨能力的情况下，可实现星上自主位置保持控制与角动量管理。

附图说明

图1为本发明实现流程图；

图2是本发明所采用的星体电推力器布局示意图；

图3是本发明中电推力器点火位置示意图；

图4是本发明中控制周期内电推力器点火时序示意图；

图5为实施例中一年内东西方向的位置保持控制结果；

图6为实施例中一年内南北方向的位置保持控制结果；

图7为实施例中一年内角动量卸载控制结果。

具体实施方式

本发明原理：步骤一可以根据测定轨所需时间以及位置保持精度要求等外部输入条件直接确定位置保持周期等参数，避免了参数选择的盲目性；步骤二根据测轨数据计算得到轨道要素控制量，其计算过程需要太阳平黄经、相对儒略日和过春分点时刻等参数，需要地面或星上提供，平经度控制量的计算过程采用了坐标轮换搜索算法，求解过程需要迭代；步骤三是计算各电推力器点火参数：点火起始时刻和点火时长，计算过程需要步骤二解得的轨道要素控制量作为输入，以及电推力器推力、推力投影系数和航天器重量等参数；步骤四是根据角动量卸载需求计算得到电推力器转动角度，计算过程需要步骤三解得的电推力器点火参数作为输入，其实质是求解一个二次规划问题，计算过程需要对初次计算结果的二阶误差进行补偿，因此实际要进行两次求解。步骤四需要每个控制周期执行一次，每个位保周期内需执行m次，步骤二和步骤三在每个位保周期内只需执行一次。执行过程如图1所示，图中实线框流程表示本发明内容对应的实际控制过程，虚线框流程表示与本发明相关的实际控制过程，实线循环表示位保周期循环，点划线循环表示控制周期循环。

下面结合具体的实例对本发明做进一步说明。

本发明的同步轨道电推进位置保持与角动量卸载联合控制方法，主要包括：确定位置保持参数、轨道要素控制量确定、电推力器点火参数确定和角动量卸载参数确定等步骤。卫星电推力器布局如图2所示，四台电推力器安装在卫星背地板，按安装位置区分为NW(1)、NE(2)、SW(3)、SE(4)，电推力器初始安装角度为θ＝50°，α＝20°，电推力器沿轨道切向、径向、法向的投影系数分别为：K_T＝0.2620、K_R＝0.7198、K_N＝0.6428。卫星初始重量3000kg，有效面积75m²，光压系数1.5。

(1)步骤一、确定位置保持周期T_s、控制周期T_c和测轨时间T_m。南北方向位置保持精度要求α_NS＝0.05°，东西方向位置保持精度要求α_EW＝0.05°。测定轨所需时间为36小时，故测轨时间T_m定为2个轨道周期即48小时，控制周期T_c定为两个轨道周期。偏心率漂移量Δe_pert可确定为：

${\mathrm{\Δe}}_{pert}=max\left(T_m,\frac{1}{2}day\right)n_s{R}_e=1.42\×{10}^{-5}$

偏心率测量误差为Δe_meas＝6×10^-6，偏心率控制误差Δe_ctrl＝1×10^-5，不确定误差Δe_rand＝0.9×10^-4，故偏心率控制精度为：

e_d＝Δe_meas+Δe_ctrl+Δe_rand+Δe_pert＝1.2×10^-4＝0.0069°

平经度控制精度可表示为：

${\mathrm{\λ}}_d={\mathrm{\λ}}_{meas}+{\mathrm{\λ}}_{ctrl}+\frac{3{\mathrm{\πK}}_T}{4K_N}{\mathrm{m\Δi}}_C+\frac{{\stackrel{\·\·}{\mathrm{\λ}}}_s{T}_m{T}_s}{16}$

位保周期需满足东西方向位保精度要求：

λ_meas、λ_ctrl分别为0.002°、0.001°，上式对m求根解得m＝6.087，向下取整得m＝6，即位保周期T_s＝T_m+mT_c＝14天，包含6个控制周期。

(2)步骤二、位置保持周期开始之前，得到一组测轨数据：(i_x0＝0，i_y0)＝(0,0.01°)，偏心率矢量(e_x0，e_y0)＝(0,0)、平经度λ₀＝-0.021°、平经度漂移率D₀＝1.81×10^-10，由倾角控制律得到每个控制周期的倾角控制量：

Δi_xC＝0.0002°

Δi_yC＝-0.0067°

由偏心率控制律得到每个控制周期的偏心率控制量为：

Δe_x＝-0.2769×10^-5

Δe_y＝-0.0242×10^-5

由平经度控制律得到平经度漂移率控制量为：

ΔD₁＝-8.4591×10^-11

ΔD₂＝-1.7534×10^-9

其中ΔD₁是前三个控制周期的平经度漂移率控制量、ΔD₂是后三个控制周期的平经度漂移率控制量。

(3)步骤三、根据各轨道要素所需控制量确定每个控制周期的电推力器点火时长及点火置参数。点火位置如图3所示，图中圆周表示一个轨道周期，X轴表示春分点方向，Δt₁、Δt₂表示NW(1)、NE(2)两台电推力器的点火弧段，均位于约赤经90°附近，Δt₃、Δt₄表示SW(3)、SE(4)两台电推力器的点火弧段，均位于约赤经270°附近。步骤二已经求得位保周期内每个控制周期需要的轨道要素控制量(Δi_x，Δi_y)、(Δe_x，Δe_y)、ΔD₁、ΔD₂，首先确定各电推力器点火弧段中点赤经：

l₁＝91.57°

l₂＝91.57°

l₃＝271.57°

l₄＝271.57°

测轨数据历元时刻为2020年1月1日00:00:00，过春分点时刻为2020年1月1日09:17:34，相对历元秒长为35454s，点火时序如图4所示，图a表示控制周期为1个轨道周期，图b表示控制周期为2个轨道周期，图中虚线刻度表示赤经，实现刻度表示相应电推力器的点火弧段。根据图4点火时序可以将电推力器点火弧段中点赤经转换为相对历元秒长：

t₁＝55371s

t₂＝141535s

t₃＝98453s

t₄＝184853s

求得前三个控制周期各电推力器所需产生的速度增量以及相应的点火时长：

$\begin{array}{c}{\mathrm{\ΔV}}_1=0.1347m/s\\ {\mathrm{\ΔV}}_2=0.1424m/s\\ {\mathrm{\ΔV}}_3=0.1438m/s\\ {\mathrm{\ΔV}}_4=0.1346m/s\end{array}\⇒\begin{array}{c}{\mathrm{\Δt}}_1=5081.0s\\ {\mathrm{\Δt}}_2=5076.1s\\ {\mathrm{\Δt}}_3=5375.1s\\ {\mathrm{\Δt}}_4=5426.9s\end{array}$

后三个控制周期各电推力器所需产生的速度增量以及相应的点火时长：

$\begin{array}{c}{\mathrm{\ΔV}}_1=0.0110m/s\\ {\mathrm{\ΔV}}_2=0.1672m/s\\ {\mathrm{\ΔV}}_3=0.1190m/s\\ {\mathrm{\ΔV}}_4=0.1593m/s\end{array}\⇒\begin{array}{c}{\mathrm{\Δt}}_1=4138.9s\\ {\mathrm{\Δt}}_2=6022.6s\\ {\mathrm{\Δt}}_3=6327.0s\\ {\mathrm{\Δt}}_4=4481.6s\end{array}$

再求得电推力器点火起始时刻(以第1个控制周期为例)：

t_start1＝52830s

t_start2＝138997s

t_start3＝95765s

t_start4＝181903s

(4)步骤四、每个控制周期开始之前，根据角动量所需卸载量确定电推力器指向所需转动的角度。以第一个控制周期的角动量卸载过程求解为例：所需卸载角动量ΔH_C＝(5,5,5)^TNms，求解模型为：

minJ＝X^TX

s.t.ΔH_C＝AX

X＝(Δα₁,Δα₂,Δα₃,Δα₄,Δθ₁,Δθ₂,Δθ₃,Δθ₄)^T

电推力器作用点到质心的距离为2m，其中矩阵A可得：

$A=\left[\begin{array}{cccccccc}-766.86& -794.69& 751.00& 818.39& -125.84& 146.18& 154.69& -134.30\\ 255.56& 296.87& -314.14& -272.73& -377.61& -369.80& -391.31& -402.98\\ 0& 0& 0& 0& 474.34& 473.89& 501.46& 506.22\end{array}\right]$

采用二次规划求解算法解得上述问题的最优解为：

X＝(0.339°,0.520°,-0.550°,-0.362°,0.126°,0.158°,0.167°,0.134°)^T

计算上述解的二阶误差为：

ΔH_2d＝(0.0046,-0.0084,0.003)^TNms

得到新的补偿二阶误差的求解模型：

minJ＝X^TX

s.t.ΔH_C+ΔH_2d＝AX

X＝(Δα₁,Δα₂,Δα₃,Δα₄,Δθ₁,Δθ₂,Δθ₃,Δθ₄)^T

求得其新的最优解为：

X＝(0.339°,0.520°,-0.550°,-0.362°,0.126°，0.158°,0.167°,0.134°)^T

与之类似可得到其他控制周期所需的电推力器转角。

重复步骤二到步骤四，在电推力器控制下卫星将保持在定点位置附近东西、南北±0.05°的范围内，且每个控制周期内角动量能够持续进行卸载，不存在角动量长期积累。

图5给出了一年内东西方向的位置保持控制结果，图中横坐标为相对初始历元时刻的天数，纵坐标为卫星地理经度。

图6给出了一年内南北方向的位置保持控制结果，图中横坐标为相对初始历元时刻的天数，纵坐标为卫星地理纬度。

图7给出了一年内角动量卸载控制结果，图中横坐标为相对初始历元时刻的位保周期数，纵坐标为该位保周期内每个控制周期结束时卫星角动量(每个控制周期开始之前卫星待卸载角动量ΔH_C＝(5,5,5)^TNms)。

Claims (10)

1.一种同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于实现步骤如下：

步骤一、确定位置保持周期T_s、控制周期T_c和测轨时间T_m

控制周期T_c为包含所有推力器各完成一次点火过程的总时间；测轨时间T_m由测定轨所需时间向上取整为整数个轨道周期得到；一个位置保持周期T_s应包含一个测轨过程T_m和若干个控制周期T_c：

T_s＝T_m+mT_c

其中m为正整数；

步骤二、每个位置保持周期T_s开始之前，根据测轨得到的轨道数据计算每个控制周期轨道要素控制量：倾角控制量矢量(Δi_x,Δi_y)、偏心率控制量矢量(Δe_x,Δe_y)和平经度漂移率控制量ΔD；

步骤三、根据各轨道要素所需控制量确定每个控制周期的电推力器点火参数，点火参数包括点火时长及点火位置；

步骤四、每个控制周期T_c开始之前，根据电推力器点火参数及角动量卸载量确定电推力器转动角度。

2.根据权利要求1所述的同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于：所述步骤一中m由东西方向位置保持精度要求确定，东西方向位置保持精度α_EW由偏心率控制精度e_d和平经度控制精度λ_d决定，如下：

α_EW＝2e_d(T_m)+λ_d(T_m,m)

m由上式求根并向下取整得到。

3.根据权利要求2所述的同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于：所述偏心率控制精度e_d按下式计算：

e_d＝Δe_meas+Δe_ctrl+Δe_rand+Δe_pert

其中Δe_meas、Δe_ctrl、Δe_rand、Δe_pert分别为偏心率测量误差、偏心率控制误差、不确定误差、偏心率漂移量；不确定误差Δe_rand包含了控制律中未予考虑的其他摄动项的影响，主要是日周期项，该项约为0.9×10^-4；偏心率漂移量Δe_pert表示无控弧段偏心率自由摄动所能达到的最大值，该项约为：

${\mathrm{\Δe}}_{pert}=max(T_m,\frac{1}{2}day)n_s{R}_e$

其中n_s为太阳视运动角速率，R_e＝0.011C_rS/m为偏心率摄动圆半径，C_r为光压系数，S/m为卫星面积质量比。

4.根据权利要求2所述的同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于：所述平经度控制精度λ_d按下式计算：

${\mathrm{\λ}}_d={\mathrm{\λ}}_{meas}+{\mathrm{\λ}}_{ctrl}+\frac{3{\mathrm{\πK}}_T}{4K_N}{\mathrm{m\Δi}}_C+\frac{{\stackrel{\·\·}{\mathrm{\λ}}}_s{T}_m{T}_s}{16}$

其中λ_meas、λ_ctrl分别为平经度测量误差、平经度控制误差；Δi_C为一个控制周期内的倾角控制量，当控制周期T_c取为1天，Δi_C为0.0026°，当控制周期T_c取为2天，Δi_C为0.0052°；ω_e为地球自转角速率，为定点位置处平经度漂移加速度；K_T、K_N分别为电推力器推力在轨道切向、法向的投影系数，取正值。

5.根据权利要求1所述的同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于：所述步骤二中的倾角控制量(Δi_x,Δi_y)按照如下控制律计算：

其中为太阳平黄经，t为距春分时间；A_x＝0.0235°，A_y＝0.0215°分别为倾角矢量半年周期摄动项的幅值；Ω_ms为白道升交点黄经，其变化记为Ω_ms＝12.11°-0.052954T，其中T为相对儒略日，相对于1950年1月1日。

6.根据权利要求1所述的同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于：所述步骤二中偏心率控制量Δe_x,Δe_y按照如下控制律计算：

其中i_s＝23.45°为黄赤交角，为太阳平黄经，t为距春分点时间；R_e＝0.011C_rS/m为偏心率摄动圆半径，C_r为光压系数，S/m为卫星面积质量比；e_x0，e_y0为偏心率矢量测量值。

7.根据权利要求1所述的同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于：所述步骤二中和平经度漂移率控制量ΔD，一个位保周期内m个控制周期分成两段，前m/2个控制周期产生相同的平经度漂移率控制量ΔD₁，后m/2个控制周期产生相同的平经度漂移率控制量ΔD₂，按照如下步骤计算：

(1)计算平经度控制目标，平经度控制目标为在当前位保周期结束时平经度和平经度漂移率控制到指定值λ^*、D^*，λ^*、D^*由下式计算：

${\mathrm{\λ}}^{*}=-\frac{5}{4}{\stackrel{\·\·}{\mathrm{\λ}}}_s\·{\left(\frac{T_m}{2}\right)}^2$

$D^{*}=-{\stackrel{\·\·}{\mathrm{\λ}}}_s\·\frac{T_m}{2}$

其中为定点位置处平经度漂移加速度；

(2)求解如下无约束优化问题：

$\min J=\min \left|\right|\left[\begin{array}{cc}\left(\frac{3m^2}{8}+\frac{m}{4}\right)T_c& \left(\frac{m^2}{8}+\frac{m}{4}\right)T_c\\ \frac{m}{2}{T}_c& \frac{m}{2}{T}_c\end{array}\right]X+\left[\begin{array}{c}{\mathrm{\λ}}_0+D_0{\mathrm{mT}}_c+\frac{1}{2}{\stackrel{\·\·}{\mathrm{\λ}}}_s{\left({\mathrm{mT}}_c\right)}^2-{\mathrm{\λ}}^{*}\\ D_0{T}_c+{\stackrel{\·\·}{\mathrm{\λ}}}_s{\mathrm{mT}}_c^2-D^{*}{T}_c\end{array}\right]|{|}_2$

X＝(ΔD₁,ΔD₂)^T

其中λ₀、D₀分别为平经度、平经度漂移率的测量值，求得一组最优解X^*＝(ΔD₁,ΔD₂)^T作为平经度漂移率的控制量；

8.根据权利要求1所述的同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于：所述步骤三中电推力器点火参数计算过程如下：

(1)首先根据倾角控制量(Δi_x,Δi_y)计算各电推力器点火弧段中点赤经：

l₁＝atan2(-Δi_y,-Δi_x)

l₂＝atan2(-Δi_y,-Δi_x)

l₃＝atan2(-Δi_y,-Δi_x)+π

l₄＝atan2(-Δi_y,-Δi_x)+π

(2)计算各电推力器所需产生的速度增量：

${\mathrm{\ΔV}}_1=\frac{1}{4}(\frac{V_s\mathrm{\Δ}i}{K_N}+\frac{V_s\left({\mathrm{\Δe}}_x\cos l_1+{\mathrm{\Δe}}_y\sin l_1\right)}{2K_T}+\frac{V_s(-{\mathrm{\Δe}}_x\sin l_1+{\mathrm{\Δe}}_y\cos l_1)}{K_R}-\frac{R_s\mathrm{\Δ}D}{3K_T})$

${\mathrm{\ΔV}}_2=\frac{1}{4}(\frac{V_s\mathrm{\Δ}i}{K_N}-\frac{V_s\left({\mathrm{\Δe}}_x\cos l_1+{\mathrm{\Δe}}_y\sin l_1\right)}{2K_T}+\frac{V_s(-{\mathrm{\Δe}}_x\sin l_1+{\mathrm{\Δe}}_y\cos l_1)}{K_R}+\frac{R_s\mathrm{\Δ}D}{3K_T})$

${\mathrm{\ΔV}}_3=\frac{1}{4}(\frac{V_s\mathrm{\Δ}i}{K_N}+\frac{V_s\left({\mathrm{\Δe}}_x\cos l_1+{\mathrm{\Δe}}_y\sin l_1\right)}{2K_T}+\frac{V_s(-{\mathrm{\Δe}}_x\sin l_1+{\mathrm{\Δe}}_y\cos l_1)}{K_R}-\frac{R_s\mathrm{\Δ}D}{3K_T})$

${\mathrm{\ΔV}}_4=\frac{1}{4}(\frac{V_s\mathrm{\Δ}i}{K_N}+\frac{V_s\left({\mathrm{\Δe}}_x\cos l_1+{\mathrm{\Δe}}_y\sin l_1\right)}{2K_T}-\frac{V_s(-{\mathrm{\Δe}}_x\sin l_1+{\mathrm{\Δe}}_y\cos l_1)}{K_R}+\frac{R_s\mathrm{\Δ}D}{3K_T})$

其中K_T、K_R、K_N分别为卫星对地姿态下电推力器推力沿轨道切向、径向、法向的投影系数，均取正值；

(3)根据各电推力器所需产生的速度增量计算实际的电推力器点火时长：

${\mathrm{\Δt}}_i=\frac{2}{{\mathrm{\ω}}_e}\arcsin \left(\frac{m_c{\mathrm{\ω}}_e{\mathrm{\ΔV}}_i}{2F_p}\right)$

其中m_c为航天器重量，F_p为电推力器推力；

(4)根据电推力器点火时长与电推力器点火弧段中点赤经计算电推力器起始点火时刻：

$t_{start}=t_0+\frac{l_i}{{\mathrm{\ω}}_e}-\frac{1}{2}{\mathrm{\Δt}}_i$

其中t₀为过春分点时刻，ω_e为地球自转角速率。

9.根据权利要求1所述的同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于：所述步骤四中根据电推力器点火参数及角动量卸载量确定电推力器转动角度，计算过程如下：

(1)求解如下二次规划问题：

minJ＝X^TX

s.t.ΔH_C＝AX

X＝(Δθ₁,Δθ₂,Δθ₃,Δθ₄,Δα₁,Δα₂,Δα₃,Δα₄)^T

其中A∈R^3×8，A的第k列为：

$A_k=\left\{\begin{array}{c}{m}_{c}d\underset{i=1}{\overset{4}{\Σ}}{\mathrm{\ΔV}}_i{C}_x\left(\frac{\mathrm{\π}}{2}\right)C_y\left(l_i+\frac{\mathrm{\π}}{2}\right)\left[\begin{array}{c}-{\mathrm{cos\α}}_i\\ 0\\ {\mathrm{sin\α}}_i\end{array}\right],k=1,2,3,4\\ m_{c}d\underset{i=1}{\overset{4}{\Σ}}{\mathrm{\ΔV}}_i{C}_x\left(\frac{\mathrm{\π}}{2}\right)C_y\left(l_i+\frac{\mathrm{\π}}{2}\right)\left[\begin{array}{c}{\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i{\mathrm{sin\α}}_i\\ {\sin }^2{\mathrm{\θ}}_i\\ {\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i{\mathrm{cos\α}}_i\end{array}\right],k=5,6,7,8\end{array}\right.$

其中θ为推力方向与Y轴的夹角；α为推力在XOZ面的投影与Z轴的夹角，α_i、θ_i分别为各电推力器初始安装角度；d为电推力器作用点到质心的距离，C_x为绕坐标系X轴旋转的方向余弦矩阵，C_y为绕坐标系Y轴旋转的方向余弦矩阵，α_i、θ_i分别为各电推力器初始安装角度。解得一组最优解X⁰＝(Δθ₁,Δθ₂,Δθ₃,Δθ₄,Δα₁,Δα₂,Δα₃,Δα₄)^T为电推力器转角的初值；

(2)根据电推力器转角初值X⁰计算角动量卸载截断误差ΔH_2d：

${\mathrm{\ΔH}}_{2d}=m_c{\mathrm{d\Σ\ΔV}}_i{C}_x\left(\frac{\mathrm{\π}}{2}\right)C_y\left(l_i+\frac{\mathrm{\π}}{2}\right)\left(-\frac{1}{2}\left[\begin{array}{c}-{\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i{\mathrm{cos\α}}_i\\ 0\\ {\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i{\mathrm{sin\α}}_i\end{array}\right]{\mathrm{\Δ\α}}_i^2+\left[\begin{array}{c}{\cos }^2{\mathrm{\θ}}_i{\mathrm{sin\α}}_i\\ -{\mathrm{cos\θ}}_i{\mathrm{sin\θ}}_i\\ {\cos }^2{\mathrm{\θ}}_i{\mathrm{cos\α}}_i\end{array}\right]{\mathrm{\Δ\θ}}_i{\mathrm{\Δ\α}}_i\right)$

(3)考虑角动量卸载截断误差ΔH_2d重新求解如下二次规划问题：

minJ＝X^TX

s.t.ΔH_C+ΔH_2d＝AX

X＝(Δα₁,Δα₂,Δα₃,Δα₄,Δθ₁,Δθ₂,Δθ₃,Δθ₄)^T

得到一组最优解X^*＝(Δθ₁,Δθ₂,Δθ₃,Δθ₄,Δα₁,Δα₂,Δα₃,Δα₄)^T作为最终的电推力器转角值。

10.根据权利要求1所述的同步轨道电推进位置保持与角动量卸载联合控制方法，其特征在于：所述步骤一中控制周期T_c可取1天或2天。