CN103991559A - 一种洛伦兹航天器悬停控制方法 - Google Patents

Abstract

传统航天器轨道控制一般采用化学推进方式，因而任务周期受所携推进工质质量制约，尤其是对能量需求较高的低轨悬停控制，所需能量之多限制了采用传统推进方式对低轨空间目标长期悬停的可实现性。本发明提出了一种采用带电洛伦兹航天器在地磁场中运动时受到的洛伦兹力与传统推力器推力联合控制的燃料最优悬停控制方法，创新性地解决了这一问题，将悬停可实现范围延伸至低轨。本发明基于建立的地球非球形摄动作用下的洛伦兹航天器相对运动模型，给出了燃料最优条件下实现悬停所需的洛伦兹航天器最优荷质比与推力器提供的最优控制加速度。

Description

一种洛伦兹航天器悬停控制方法

技术领域

本发明涉及航天器悬停控制方法，特别涉及一种采用洛伦兹力与传统推力联合进行悬停控制的方法。

背景技术

所谓航天器悬停，是指通过力和力矩控制，保持航天器相对于某空间目标的位置始终保持不变，航天器相对于空间目标仿佛是静止“悬停”于某个固定点上。悬停轨道在空间任务中的应用广泛。若将航天器悬停视为一种特殊形式的空间编队飞行，则最简单的悬停构型为跟飞构型，该构型中，追踪器与目标器位于同一圆轨道，仅提前或落后目标一固定相位角。除该特殊悬停构型，为实现悬停，需对航天器持续施加控制加速度以保持其对目标的相对位置不变。因而，悬停轨道多为非开普勒轨道。传统航天器轨道控制一般采用化学推进方式。因而，空间任务周期往往受航天器所携推进工质质量制约，尤其是对能量需求较高的低轨悬停控制。例如，对500km赤道圆轨道的目标航天器在其径向上方5km悬停一天所需的速度增量为1.586km/s。所需能量之多限制了采用传统推进方式对低轨空间目标长期悬停的可实现性，成为制约低轨航天器执行悬停任务的技术难点。

空间带电洛伦兹航天器在地磁场中运动时通常受到洛伦兹力作用，该洛伦兹力可用以进行航天器绝对轨道与相对轨道控制，从而实现无推进工质消耗的轨道机动。洛伦兹力由当地磁场以及航天器与当地磁场相对速度决定。由于低轨(轨道高度：100km至2000km)磁场强度明显强于高轨(轨道高度：20000km以上)且相对速度更快，因而，洛伦兹航天器在低轨应用效率更高。采用洛伦兹力辅助悬停将有效减少燃料消耗，使得长期悬停任务的可实现范围延伸至低轨。

发明内容

本发明针对洛伦兹力辅助悬停控制问题，假设地磁场为一随地球自转的倾斜磁偶极子且洛伦兹航天器可视为带电质点，采用拉格朗日动力学方法建立了地球非球形摄动二阶带谐系数(以下简称J₂)摄动作用下的洛伦兹力辅助悬停动力学模型，基于所建立的悬停动力学模型，分析给出了实现悬停所需的开环总控制加速度，由于洛伦兹力的方向始终与当地磁场方向以及航天器相对当地磁场速度方向垂直，这一特性决定了洛伦兹力方向并不一定与悬停所需开环总控制加速度方向相同。因此，需要采用推力器推力与洛伦兹力联合控制，其中洛伦兹力作为辅助推力，减少推进工质的消耗。如何确定洛伦兹力与推力器推力的最优配比从而使得悬停任务周期内实现悬停所需速度增量消耗最少，是本发明设计中的技术难点。本发明基于欧拉-拉格朗日方程提出了一种洛伦兹力与推力器推力联合控制的燃料最优混合动力悬停控制方法。

本发明的技术方案如下：

一种航天器悬停控制方法，按以下步骤进行：

步骤一：给定初始时刻目标轨道参数与地磁场参数，给定悬停构型；

步骤二：各轴悬停距离计算：计算追踪器与目标器的相对位置矢量在相对运动坐标系中的各轴向分量；

步骤三：悬停所需总控制加速度计算：推导J₂摄动作用下洛伦兹力辅助悬停动力学方程，根据悬停定义，计算实现悬停所需的各轴向总控制加速度；

步骤四：悬停所需洛伦兹航天器最优荷质比计算：以任务历时内实现悬停所需的速度增量消耗最优为指标，求解欧拉-拉格朗日方程得到该指标下洛伦兹航天器最优荷质比轨迹；

步骤五：悬停所需最优洛伦兹加速度及推力器提供的最优控制加速度计算：将洛伦兹航天器最优荷质比代入洛伦兹加速度定义式得到最优洛伦兹加速度，通过力的分解得到实现悬停所需的由推力器提供的最优控制加速度。

由于在地球静止或地球同步轨道(轨道高度约为35786km、轨道周期等于地球自转周期的圆轨道，以下简称GEOs)，洛伦兹航天器与地磁场不存在相对运动速度，因而无法产生洛伦兹力，也就无法对此类特殊轨道目标进行悬停，除去该类特殊轨道外，洛伦兹航天器均可适用；考虑到洛伦兹航天器在低轨应用效率高，且对低轨目标悬停所需速度增量大，优选的，本发明应用于对低轨空间目标的悬停，即轨道高度在100km至2000km；

在步骤一中所述的悬停构型Γ＝[ρ,E,A]，其中ρ为悬停距离，E与A分别为悬停高低角与悬停方位角，其定义方式如图1所示，图中，O_EX_IY_IZ_I为地心惯性坐标系，O_E为地心，相对运动在相对运动坐标系O_Txyz中描述，该坐标系原点位于目标器质心O_T，其中x轴沿目标器径向，z轴与目标器轨道面法向重合，y轴与x、z轴构成笛卡尔右手直角坐标系，O_L为洛伦兹航天器质心，悬停高低角E的范围为E∈[-π/2,π/2]，规定追踪器悬停于目标器上方时为正，反之为负；悬停方位角A的范围为A∈[-π,π]，从飞行方向(+y)看，规定追踪器悬停于目标器左方时为正，反之为负，根据上述定义，图中的悬停高低角与方位角均为正；

在步骤二中所述的计算追踪器与目标器的相对位置矢量在相对运动坐标系中的各轴向分量，其计算方法为：

$\left\{\begin{array}{c}x=\mathrm{\ρ}\sin E\\ y=\mathrm{\ρ}\cos E\cos A\\ z=\mathrm{\ρ}\cos E\sin A\end{array}\right.---\left(1\right)$

式中，ρ＝[x y z]^T为追踪器与目标器的相对位置矢量在相对运动坐标系中的描述，上标T为矩阵转置符号；

在步骤三中所述的实现悬停所需的各轴向总控制加速度，其计算方法如下：

(1)建立J₂摄动作用下的洛伦兹航天器相对运动动力学模型：

基于拉格朗日动力学方法推导得到的J₂摄动作用下的洛伦兹力辅助悬停动力学方程：

$\left\{\begin{array}{c}\stackrel{\·\·}{x}=2\stackrel{\·}{y}{\mathrm{\ω}}_z-x({\mathrm{\η}}_L^2-{\mathrm{\ω}}_z^2)+{\mathrm{y\ϵ}}_z-{\mathrm{z\ω}}_x{\mathrm{\ω}}_z-({\mathrm{\ξ}}_L-\mathrm{\ξ})\sin i_T{\sin u}_T-R_T({\mathrm{\η}}_L^2-{\mathrm{\η}}^2)+a_x+a_R\\ \stackrel{\·\·}{y}=-2\stackrel{\·}{x}{\mathrm{\ω}}_z+2\stackrel{\·}{z}{\mathrm{\ω}}_x-x{\mathrm{\ϵ}}_z-y({\mathrm{\η}}_L^2-{\mathrm{\ω}}_z^z-{\mathrm{\ω}}_x^2)+{\mathrm{z\ϵ}}_x-({\mathrm{\ξ}}_L-\mathrm{\ξ})\sin i_T\cos u_T+a_y+a_S\\ \stackrel{\·\·}{z}=-2\stackrel{\·}{y}{\mathrm{\ω}}_x-x{\mathrm{\ω}}_x{\mathrm{\ω}}_z-y{\mathrm{\ω}}_x-z({\mathrm{\η}}_L^2-{\mathrm{\ω}}_x^2)-({\mathrm{\ξ}}_L-\mathrm{\ξ})\cos i_T+a_z+a_W\end{array}\right.---\left(2\right)$

其中，

$\left\{\begin{array}{c}{\mathrm{\η}}_L^2=\frac{\mathrm{\μ}}{R_L^3}+\frac{k_J}{R_L^5}-\frac{{5k}_J{R}_{\mathrm{Lz}}^2}{R_L^7}\\ {\mathrm{\η}}^2=\frac{\mathrm{\μ}}{R_T^3}+\frac{k_J}{R_T^5}-\frac{{5k}_J\mathrm{si}{n}^2{i}_T{\sin }^2{u}_T}{R_T^5}\end{array}\right.---\left(3\right)$

$\left\{\begin{array}{c}{\mathrm{\ξ}}_L=\frac{{2k}_J{R}_{\mathrm{Lz}}}{R_L^5}\\ \mathrm{\ξ}=\frac{{2k}_J\sin i_T\sin u_T}{R_T^4}\end{array}\right.---\left(4\right)$

式中，洛伦兹航天器地心距为R_L＝[(R_T+x)²+y²+z²]^1/2；

ω＝[ω_x ω_y ω_z]^T与ε＝[ε_x ε_y ε_z]^T分别为相对运动坐标系的角速度与角加速度，即

$\mathrm{\ω}=\left[\begin{array}{c}{\mathrm{\ω}}_x\\ {\mathrm{\ω}}_y\\ {\mathrm{\ω}}_z\end{array}\right]=\left[\begin{array}{c}-\left(k_J\sin 2i_T\sin u_T\right)/\left(h_T{R}_T^3\right)\\ 0\\ h_T/R_T^2\end{array}\right]---\left(5\right)$

$\mathrm{\ϵ}=\left[\begin{array}{c}{\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_y\\ {\mathrm{\ϵ}}_z\end{array}\right]=\left[\begin{array}{c}\frac{{3V}_r{k}_J\sin {2i}_T\sin u_T}{R_T^4{h}_T}-\frac{k_J\sin 2i_T\cos u_T}{R_T^5}-\frac{{8k}_J^2\mathrm{si}{n}^3{i}_T\cos i_T{\sin }^2{u}_T{\cos u}_T}{R_T^6{h}_T^2}\\ 0\\ -\frac{{2h}_T{V}_r}{R_T^3}-\frac{k_J{\sin }^2{i}_T\sin {2u}_T}{R_T^5}\end{array}\right]---\left(6\right)$

式中，R_Lz＝(R_T+x)sin i_Tsin u_T+y sin i_T cos u_T+z cos i_T为洛伦兹航天器地心距在地惯系中Z_I方向的投影，i_T与u_T分别为目标器的轨道倾角与纬度幅角，系数其中地球非球形摄动二阶带谐系数为J₂＝1.0826×10^-3，μ＝3.986005×10¹⁴m³/s²为地球引力常数，R_E为地球半径；

a_C＝[a_R a_S a_W]^T与a_L＝λl＝λ[l _xl_y l_z]^T分别为作用于洛伦兹航天器的控制加速度与洛伦兹加速度，其中λ＝q/m为洛伦兹航天器的荷质比，且

$\left\{\begin{array}{c}{l}_x=\frac{B_0}{R_L^5}{V}_y[({3z}^2-R_L^2)n_z+3z[(R_T+x)n_x+{\mathrm{yn}}_y\left]\right\}-\frac{B_0}{R_L^5}{V}_{z,J_2}\{({3y}^2-R_L^2)n_y+3y[(R_T+x)n_x+{\mathrm{zn}}_z\left]\right\}\\ l_y=\frac{B_0}{R_L^5}{V}_z\left\{\right[3{(R_T+x)}^2-R_L^2]n_z+3(R_T+x)({\mathrm{yn}}_y+{\mathrm{zn}}_z)\}-\frac{B_0}{R_L^5}{V}_x[({3z}^2-R_L^2)n_z+3z[(R_T+x)n_x+{\mathrm{zn}}_z\left]\right\}\\ l_z=\frac{B_0}{R_L^5}{V}_x\{({3y}^2-R_L^2)n_y+3y[(R_T+x)n_x+{\mathrm{zn}}_x\left]\right\}-\frac{B_0}{R_L^5}{V}_y\left\{\right[3{(R_T+x)}^2-R_L^2]n_z+3(R_T+x)({\mathrm{yn}}_y+{\mathrm{zn}}_z)\}\end{array}\right.---\left(7\right)$

式中，B₀＝8.0×10¹⁵T·m³为地磁场地磁矩值。V_rel＝[V _xV_y V_z]^T为洛伦兹航天器与当地磁场相对速度在相对运动坐标系中的描述：

$V_{\mathrm{rel}}=\left[\begin{array}{c}{V}_x\\ V_y\\ V_z\end{array}\right]=\left[\begin{array}{c}{\stackrel{\·}{R}}_T+\stackrel{\·}{x}-y({\mathrm{\ω}}_z-{\mathrm{\ω}}_E\cos i_T)-{\mathrm{z\ω}}_E\sin i_T\cos u_T\\ \stackrel{\·}{y}+(R_T+x)({\mathrm{\ω}}_z-{\mathrm{\ω}}_E\cos i_T)-z({\mathrm{\ω}}_x-{\mathrm{\ω}}_E\sin i_T\sin u_T)\\ \stackrel{\·}{z}+(R_T+x){\mathrm{\ω}}_E\sin i_T\cos u_T+y({\mathrm{\ω}}_x-{\mathrm{\ω}}_E\sin i_T\sin u_T)\end{array}\right]---\left(8\right)$

$\hat{n}={\left[\begin{array}{ccc}{n}_x& n_y& n_z\end{array}\right]}^T$ 为单位磁偶极子矢量在相对运动坐标系中的描述：

$\hat{n}=\left[\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right]=\left[\begin{array}{c}-(\cos \mathrm{\β}\cos u_T+\sin \mathrm{\β}\cos i_T\sin u_T)\sin \mathrm{\α}-{\sin i}_T{\sin u}_T\cos \mathrm{\α}\\ (\cos \mathrm{\β}\sin u_T-\sin \mathrm{\β}\cos i_T\cos u_T)\sin \mathrm{\α}-{\sin i}_T{\cos u}_T\cos \mathrm{\α}\\ \sin \mathrm{\β}\sin i_T\sin \mathrm{\α}-{\cos i}_T\cos \mathrm{\α}\end{array}\right]---\left(9\right)$

式中，β＝Ω_M-Ω_T，其中Ω_M＝ω_Et+Ω₀为地磁轴在地惯系赤道平面内的相位角，ω_E为地球自转角速度，Ω₀为初始时刻地磁轴在地惯系赤道平面内的相位角，Ω_T为目标器轨道升交点赤经，倾斜地磁轴与地球自转轴的夹角为α，各角度定义及其与地惯系、相对运动坐标系位置关系如图2所示；

六个轨道参数(R_T,V_r,h_T,i_T,u_T,Ω_T)分别为目标器的轨道半径、径向速度、轨道角动量、轨道倾角、近地点幅角以及升交点赤经，其满足动力学方程

$\left\{\begin{array}{c}{\stackrel{\·}{R}}_T=V_r\\ {\stackrel{\·}{V}}_r=-\mathrm{\μ}/R_T^2+h_T^2/R_T^3-k_J(1-3{\sin }^2{i}_T{\sin }^2{u}_T)/R_T^4\\ {\stackrel{\·}{h}}_T=-\left(k_J{\sin }^2{i}_T\sin 2u_T\right)/R_T^3\\ {\stackrel{\·}{\stackrel{\·}{i}}}_T=-\left(k_J\sin {2i}_T\sin 2u_T\right)/\left({2h}_T{R}_T^3\right)\\ {\stackrel{\·}{u}}_T=h_T/R_T^2+\left({2k}_J\mathrm{co}{s}^2{i}_T\mathrm{si}{n}^2{u}_T\right)/\left(h_T{R}_T^3\right)\\ {\stackrel{\·}{\mathrm{\Ω}}}_T=-\left({2k}_J\cos i_T{\sin }^2{u}_T\right)/\left(h_T{R}_T^3\right)\end{array}\right.---\left(10\right)$

式(2)至式(10)即J₂摄动作用下的洛伦兹航天器非线性相对运动方程；

(2)计算悬停所需的开环控制总加速度：

根据悬停定义，在相对运动坐标系中，追踪器与目标器的相对位置矢量保持不变，因而相对位置各轴分量对时间的各阶导数为零，即

$\left\{\begin{array}{c}\stackrel{\·}{x}=\stackrel{\·}{y}=\stackrel{\·}{z}=0\\ \stackrel{\·\·}{x}=\stackrel{\·\·}{y}=\stackrel{\·\·}{z}=0\end{array}\right.---\left(11\right)$

将式(11)代入式(2)中得，

a_C＝h-λl              (12)

式中，h为实现悬停所需的总控制加速度，即

$h=\left[\begin{array}{c}x({\mathrm{\η}}_L^2-{\mathrm{\ω}}_z^2)-y{\mathrm{\ϵ}}_z+{\mathrm{z\ω}}_x{\mathrm{\ω}}_z+({\mathrm{\ξ}}_L-\mathrm{\ξ})\sin i_T\sin u_T+R_T({\mathrm{\η}}_L^2-{\mathrm{\η}}^2)\\ {\mathrm{x\ϵ}}_z+y({\mathrm{\η}}_L^2-{\mathrm{\ω}}_z^2-{\mathrm{\ω}}_x^2)-z{\mathrm{\ϵ}}_x+({\mathrm{\ξ}}_L-\mathrm{\ξ})\sin i_T\cos u_T\\ {\mathrm{x\ω}}_x{\mathrm{\ω}}_z+{\mathrm{y\ϵ}}_x+z({\mathrm{\η}}_L^2-{\mathrm{\ω}}_x^2)+({\mathrm{\ξ}}_L-\mathrm{\ξ})\cos i_T\end{array}\right]---\left(13\right)$

在步骤四中所述的实现悬停所需的最优洛伦兹航天器荷质比，其计算方法如下：

由于洛伦兹加速度方向与当地磁场方向及洛伦兹航天器相对磁场速度方向垂直，当洛伦兹航天器轨道位置确定，洛伦兹加速度方向也随之确定；若洛伦兹加速度方向与悬停所需加速度方向重合，即l∥h，则若洛伦兹航天器荷质比大小为λ＝||h||/||l||，产生的洛伦兹加速度即可完全抵偿悬停所需加速度；然而，由于洛伦兹力作用方向的局限性，并不能保证洛伦兹加速度方向与悬停所需加速度方向重合，因此，为实现悬停，需施加额外的控制加速度，本发明中假设该额外控制加速度由推力器提供，为最大限度地节省燃料，选取拉格朗日形式的燃料最优目标函数，即：

$J={\∫}_0^{t_f}L[t,\mathrm{\λ}\left(t\right)]\mathrm{dt}={\∫}_0^{t_f}\left|\right|a_C\left|\right|\mathrm{dt}---\left(14\right)$

式中，t_f为悬停任务历时；

求解欧拉-拉格朗日方程：

$\frac{d}{\mathrm{dt}}\left(\frac{\∂L}{\∂\stackrel{\·}{\mathrm{\λ}}}\right)-\frac{\∂L}{\∂\mathrm{\λ}}=0---\left(15\right)$

得到实现燃料最省的洛伦兹航天器荷质比最优轨迹

${\mathrm{\λ}}^{*}\left(t\right)=\left\{\begin{array}{cc}\frac{h\·l}{{\left|\right|l\left|\right|}^2},& \left|\right|l\left|\right|\≠0\\ 0,& \left|\right|l\left|\right|=0\end{array}\right.---\left(16\right)$

当磁场B方向与相对速度V_rel方向平行时，无洛伦兹加速度产生，即||l||＝0，因而，悬停所需加速度将完全由推力器提供，该条件下设定洛伦兹航天器荷质比为0，如式(16)所示；

在步骤五中所述的实现悬停所需最优洛伦兹加速度及推力器提供的最优控制加速度，其计算方法如下：

将洛伦兹航天器最优荷质比λ^*代入洛伦兹加速度定义式，得到实现悬停所需的最优洛伦兹加速度为

$a_L^{*}\left(t\right)=\left\{\begin{array}{cc}\frac{h\·l}{{\left|\right|l\left|\right|}^2}l,& \left|\right|l\left|\right|\≠0\\ 0,& \left|\right|l\left|\right|=0\end{array}\right.---\left(17\right)$

将式(17)代入式(12)，得到最优控制加速度

$a_C^{*}\left(t\right)=\left\{\begin{array}{cc}h-\frac{h\·l}{{\left|\right|l\left|\right|}^2}l,& \left|\right|l\left|\right|\≠0\\ h,& \left|\right|l\left|\right|=0\end{array}\right.---\left(18\right)$

因此，对于带电洛伦兹航天器，实现悬停所需的速度增量为

$\mathrm{\ΔV}={\∫}_0^{t_f}\left|\right|a_C^{*}\left|\right|\mathrm{dt}---\left(19\right)$

而对于非带电航天器，实现悬停所需的控制加速度均由推力器提供，因而所需速度增量为

$\mathrm{\ΔV}={\∫}_0^{t_f}\left|\right|h\left|\right|\mathrm{dt}---\left(20\right)$

与现有技术相比，本发明具有以下优点：

1)本发明直接基于洛伦兹航天器相对运动非线性动力学模型设计，考虑了航天器相对运动的非线性特性，克服了传统相对运动模型由线性化引起的误差。基此模型求解得到的悬停控制加速度更为精确。

2)本发明考虑到J₂摄动是低地球轨道的主要摄动力且悬停构型保持对摄动力较为敏感，建立了J₂摄动力与洛伦兹力共同作用下的洛伦兹航天器相对运动模型，更符合低地球轨道空间环境特性，因而求解得到的悬停控制加速度更为精确。

3)本发明采用洛伦兹力与推力器推力联合进行航天器悬停控制，创新性地解决了采用传统化学推进方式由于推进工质质量制约而无法对低轨道空间目标长期悬停的问题，利用洛伦兹力作为辅助推力，节省了速度增量消耗，延长了悬停任务周期，将长期悬停任务可实现范围延伸至低地球轨道。

控制工程师在应用过程中可以根据实际空间目标给定任意悬停构型，并将由该方法得到的控制量传输至执行机构实现对空间目标的悬停。

附图说明

图1为本发明坐标系与悬停构型定义图；

图2为本发明倾斜偶极子空间指向与相关角度定义图；

图3为本发明所述控制方法流程图；

图4为本发明实现悬停一天所需总控制加速度图；

图5为本发明实现悬停一个周期所需总控制加速度图；

图6为本发明实现悬停一天所需洛伦兹航天器最优荷质比图；

图7为本发明实现悬停一个周期所需洛伦兹航天器最优荷质比图；

图8为本发明实现悬停一天所需最优洛伦兹加速度图；

图9为本发明实现悬停一个周期所需最优洛伦兹加速度图；

图10为本发明实现悬停一天所需最优控制加速度对比图；

图11为本发明实现悬停一个周期所需最优控制加速度对比图；

图中及本发明中符号说明如下：

A悬停方位角，deg或rad

a_C控制加速度矢量，m/s²

最优控制加速度矢量，m/s²

||a_C||控制加速度大小，m/s²

a_L洛伦兹加速度矢量，m/s²

最优洛伦兹加速度矢量，m/s²

||a_L||洛伦兹加速度大小，m/s²

a_R控制加速度径向分量，m/s²

a_S控制加速度迹向分量，m/s²

a_W控制加速度法向分量，m/s²

a_x洛伦兹加速度径向分量，m/s²

a_y洛伦兹加速度迹向分量，m/s²

a_z洛伦兹加速度法向分量，m/s²

B磁场强度矢量，T

B_x磁场强度径向分量，T

B_y磁场强度迹向分量，T

B_z磁场强度法向分量，T

B₀地磁场地磁矩值，T·m³

E悬停高低角，deg或rad

h悬停所需总加速度矢量，m/s²

||h||悬停所需总加速度大小，m/s²

h_T目标航天器轨道角动量，m²/s

h_x悬停所需总加速度径向分量，m/s²

h_y悬停所需总加速度迹向分量，m/s²

h_z悬停所需总加速度法向分量，m/s²

i_T目标航天器轨道倾角，deg或rad

J指标函数，m/s

J₂地球非球形摄动二阶带谐系数

L拉格朗日函数

l洛伦兹加速度方向矢量(非单位矢量)，(m·T)/s

l_x洛伦兹加速度方向矢量径向分量，(m·T)/s

l_y洛伦兹加速度方向矢量迹向分量，(m·T)/s

l_z洛伦兹加速度方向矢量法向分量，(m·T)/s

m洛伦兹航天器质量，kg

单位偶极子矢量

n_x单位偶极子矢量径向分量

n_y单位偶极子矢量迹向分量

n_z单位偶极子矢量法向分量

O_EX_IY_IZ_I地心惯性坐标系

O_L洛伦兹航天器质心

O_T目标航天器质心

O_Txyz相对运动坐标系

q洛伦兹航天器带电量，C

q/m洛伦兹航天器荷质比，C/kg

R_E地球半径，m

R_L洛伦兹航天器地心距，m

R_Lz洛伦兹航天器地心距在地心惯性坐标系O_EX_IY_IZ_I中O_EZ_I轴上的投影，m

R_T目标航天器地心距，m

t时间，s

t_f悬停任务历时，s

u_T目标航天器纬度幅角，deg或rad

V_r目标航天器径向速度，m/s

V_rel洛伦兹航天器与当地磁场相对速度矢量，m/s

V_x洛伦兹航天器与当地磁场相对速度径向分量，m/s

V_y洛伦兹航天器与当地磁场相对速度迹向分量，m/s

V_z洛伦兹航天器与当地磁场相对速度法向分量，m/s

x追踪器与目标器径向相对距离，m

y追踪器与目标器迹向相对距离，m

z追踪器与目标器法向相对距离，m

α地磁轴倾角，deg或rad

β地惯系赤道面内地磁轴相位角与目标器轨道升交点赤经之差，deg或rad

ΔV实现悬停所需速度增量，m/s

ε相对运动坐标系转动角加速度矢量，rad/s²

ε_x相对运动坐标系转动角加速度矢量径向分量，rad/s²

ε_y相对运动坐标系转动角加速度矢量迹向分量，rad/s²

ε_z相对运动坐标系转动角加速度矢量法向分量，rad/s²

η系数

η_L系数

Γ悬停构型

λ洛伦兹航天器荷质比，C/kg

λ^*最优洛伦兹航天器荷质比，C/kg

μ地球引力常数，m³/s²

ρ追踪器与目标器相对位置矢量，m

ρ追踪器与目标器相对距离，m

Ω_M地磁轴在地惯系赤道面内相位角，deg或rad

Ω_T目标航天器轨道升交点赤经，deg或rad

Ω₀初始时刻地磁轴在地惯系赤道面内相位角，deg或rad

ω相对运动坐标系转动角速度矢量，rad/s

ω_E地球自转角速度，rad/s

ω_x相对运动坐标系转动角速度矢量径向分量，rad/s

ω_y相对运动坐标系转动角速度矢量迹向分量，rad/s

ω_z相对运动坐标系转动角速度矢量法向分量，rad/s

ξ系数

ξ_L系数

具体实施方案

下面结合附图，对本发明中的设计方法作进一步的说明：

本发明“一种航天器悬停控制方法”，其具体步骤如下：

步骤一：给定初始时刻目标轨道参数、地磁场参数与悬停构型：

如图3中所示，给定初始时刻目标航天器轨道参数如表1所示，目标轨道周期约为1.6h，属于低轨近圆轨道，初始时刻地磁轴在地惯系赤道面内的相位角Ω₀＝-60°，悬停构型为Γ＝[4km,π/2,0]，其中悬停距离ρ为4km，悬停高低角E为90°，悬停方位角A为0°；

表1初始时刻目标航天器轨道根数

步骤二：各轴悬停距离计算

如图3中所示，计算悬停距离在相对运动坐标系各轴向分量。将悬停构型代入式(1)，得到各轴向的悬停距离分别为：

$\left\{\begin{array}{c}x=4\mathrm{km}\\ y=0\mathrm{km}\\ z=0\mathrm{km}\end{array}\right.---\left(21\right)$

步骤三：悬停所需总控制加速度计算

如图3中所示，计算实现悬停所需的各轴向总控制加速度。将步骤一与步骤二中的参数代入式(13)，得到悬停任务过程中任一时刻所需的总控制加速度h，其中，目标航天器的六个轨道参数(R_T,V_r,h_T,i_T,u_T,Ω_T)为对式(10)进行数值积分的结果，一天内实现该悬停构型所需的各轴向总控制加速度如图4所示，一个轨道轨道周期内实现该悬停构型所需的各轴向总控制加速度如图5所示；

步骤四：悬停所需洛伦兹航天器最优荷质比计算

如图3中所示，计算燃料最优指标条件下实现悬停所需的洛伦兹航天器最优荷质比。将步骤一与步骤二中的参数代入式(7)中，得到悬停任务过程总任一时刻洛伦兹加速度的方向矢量。将该方向矢量l与悬停所需总加速度矢量h代入式(16)，即得到悬停任务过程中任一时刻实现悬停所需的洛伦兹航天器最优荷质比，如图6与图7所示所示，其中，图6为一天内的最优洛伦兹航天器荷质比，图7为一个周期内的最优洛伦兹航天器荷质比；

步骤五：悬停所需最优洛伦兹加速度及推力器提供的最优控制加速度计算

如图3中所示，计算燃料最优指标条件下实现悬停所需的最优洛伦兹加速度以及由推力器提供的最优控制加速度；将步骤四中得到的洛伦兹航天器最优荷质比λ^*与洛伦兹加速度方向矢量l代入式(17)，即得到实现悬停所需的最优洛伦兹加速度，如图8和图9所示，其中，图8为一天内的最优洛伦兹加速度，图9为一个周期内的最优洛伦兹加速度，同理，将步骤三中得到的悬停所需总加速度矢量h与步骤四中得到的洛伦兹加速度方向矢量l代入式(18)，即得到采用洛伦兹航天器实现悬停所需的由推力器提供的最优控制加速度如图10和图11所示，图10中给出了一天内采用传统非带电航天器与带电洛伦兹航天器实现悬停所需的由推力器提供的控制加速度对比图，同理，图11为一个周期内的对比图。

将实现悬停所需的总控制加速度h轨迹代入式(20)进行数值积分，得到非带电航天器悬停所需的速度增量，悬停一个周期所需的增量为82.71m/s，悬停一天(按24h计)所需的速度增量为1.24×10³m/s，将洛伦兹航天器实现悬停所需的控制加速度轨迹代入式(19)进行数值积分，得到采用洛伦兹力辅助悬停所需的速度增量，采用洛伦兹力辅助悬停，一个周期内所需速度增量为4.54m/s，节省了约94.51％速度增量，一天内所需速度增量为365.85m/s，节省了约70.51％速度增量，对比发现，采用洛伦兹力辅助悬停，可有效减少速度增量消耗，延长悬停任务历时，将长期悬停可实现范围延伸至低地球轨道。

Claims (3)

1.一种航天器悬停控制方法，应用于非地球静止或地球同步轨道空间目标的悬停，按以下步骤进行：

步骤一：给定初始时刻目标轨道参数与地磁场参数，给定悬停构型；

步骤二：各轴悬停距离计算：计算追踪器与目标器的相对位置矢量在相对运动坐标系中的各轴向分量；

步骤三：悬停所需总控制加速度计算：推导J₂摄动作用下洛伦兹力辅助悬停动力学方程，根据悬停定义，计算实现悬停所需的各轴向总控制加速度；

步骤四：悬停所需洛伦兹航天器最优荷质比计算：以任务历时内实现悬停所需的速度增量消耗最优为指标，求解欧拉-拉格朗日方程得到该指标下洛伦兹航天器最优荷质比轨迹；

步骤五：悬停所需最优洛伦兹加速度及推力器提供的最优控制加速度计算：将洛伦兹航天器最优荷质比代入洛伦兹加速度定义式得到最优洛伦兹加速度，通过力的分解得到实现悬停所需的由推力器提供的最优控制加速度。

2.权利要求1所述一种航天器悬停控制方法，其特征在于：该方法应用于对低轨空间目标的悬停，即轨道高度在100km至2000km。

3.权利要求1或2所述一种航天器悬停控制方法，其特征在于：

在步骤一中所述的悬停构型Γ＝[ρ,E,A]，其中ρ为悬停距离，E与A分别为悬停高低角与悬停方位角，其定义方式如图1所示，O_EX_IY_IZ_I为地心惯性坐标系，O_E为地心，相对运动在相对运动坐标系O_Txyz中描述，该坐标系原点位于目标器质心O_T，其中x轴沿目标器径向，z轴与目标器轨道面法向重合，y轴与x、z轴构成笛卡尔右手直角坐标系，O_L为洛伦兹航天器质心，悬停高低角E的范围为E∈[-π/2,π/2]，规定追踪器悬停于目标器上方时为正，反之为负；悬停方位角A的范围为A∈[-π,π]，从飞行方向(+y)看，规定追踪器悬停于目标器左方时为正，反之为负，根据上述定义，图中的悬停高低角与方位角均为正；

在步骤二中所述的计算追踪器与目标器的相对位置矢量在相对运动坐标系中的各轴向分量，其计算方法为：

$\left\{\begin{array}{c}x=\mathrm{\ρ}\sin E\\ y=\mathrm{\ρ}\cos E\cos A\\ z=\mathrm{\ρ}\cos E\sin A\end{array}\right.---\left(1\right)$

式中，ρ＝[x y z]^T为追踪器与目标器的相对位置矢量在相对运动坐标系中的描述，上标T为矩阵转置符号；

在步骤三中所述的实现悬停所需的各轴向总控制加速度，其计算方法如下：

(1)建立J₂摄动作用下的洛伦兹航天器相对运动动力学模型：

基于拉格朗日动力学方法推导得到的J₂摄动作用下的洛伦兹力辅助悬停动力学方程：

$\left\{\begin{array}{c}\stackrel{\·\·}{x}=2\stackrel{\·}{y}{\mathrm{\ω}}_z-x({\mathrm{\η}}_L^2-{\mathrm{\ω}}_z^2)+{\mathrm{y\ϵ}}_z-{\mathrm{z\ω}}_x{\mathrm{\ω}}_z-({\mathrm{\ξ}}_L-\mathrm{\ξ})\sin i_T{\sin u}_T-R_T({\mathrm{\η}}_L^2-{\mathrm{\η}}^2)+a_x+a_R\\ \stackrel{\·\·}{y}=-2\stackrel{\·}{x}{\mathrm{\ω}}_z+2\stackrel{\·}{z}{\mathrm{\ω}}_x-x{\mathrm{\ϵ}}_z-y({\mathrm{\η}}_L^2-{\mathrm{\ω}}_z^z-{\mathrm{\ω}}_x^2)+{\mathrm{z\ϵ}}_x-({\mathrm{\ξ}}_L-\mathrm{\ξ})\sin i_T\cos u_T+a_y+a_S\\ \stackrel{\·\·}{z}=-2\stackrel{\·}{y}{\mathrm{\ω}}_x-x{\mathrm{\ω}}_x{\mathrm{\ω}}_z-y{\mathrm{\ω}}_x-z({\mathrm{\η}}_L^2-{\mathrm{\ω}}_x^2)-({\mathrm{\ξ}}_L-\mathrm{\ξ})\cos i_T+a_z+a_W\end{array}\right.---\left(2\right)$

其中，

$\left\{\begin{array}{c}{\mathrm{\η}}_L^2=\frac{\mathrm{\μ}}{R_L^3}+\frac{k_J}{R_L^5}-\frac{{5k}_J{R}_{\mathrm{Lz}}^2}{R_L^7}\\ {\mathrm{\η}}^2=\frac{\mathrm{\μ}}{R_T^3}+\frac{k_J}{R_T^5}-\frac{{5k}_J\mathrm{si}{n}^2{i}_T{\sin }^2{u}_T}{R_T^5}\end{array}\right.---\left(3\right)$

$\left\{\begin{array}{c}{\mathrm{\ξ}}_L=\frac{{2k}_J{R}_{\mathrm{Lz}}}{R_L^5}\\ \mathrm{\ξ}=\frac{{2k}_J\sin i_T\sin u_T}{R_T^4}\end{array}\right.---\left(4\right)$

式中，洛伦兹航天器地心距为R_L＝[(R_T+x)²+_y ²+z²]^1/2；

ω＝[ω_x ω_y ω_z]^T与ε＝[ε_x ε_y ε_z]^T分别为相对运动坐标系的角速度与角加速度，即

$\mathrm{\ω}=\left[\begin{array}{c}{\mathrm{\ω}}_x\\ {\mathrm{\ω}}_y\\ {\mathrm{\ω}}_z\end{array}\right]=\left[\begin{array}{c}-\left(k_J\sin 2i_T\sin u_T\right)/\left(h_T{R}_T^3\right)\\ 0\\ h_T/R_T^2\end{array}\right]---\left(5\right)$

$\mathrm{\ϵ}=\left[\begin{array}{c}{\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_y\\ {\mathrm{\ϵ}}_z\end{array}\right]=\left[\begin{array}{c}\frac{{3V}_r{k}_J\sin {2i}_T\sin u_T}{R_T^4{h}_T}-\frac{k_J\sin 2i_T\cos u_T}{R_T^5}-\frac{{8k}_J^2\mathrm{si}{n}^3{i}_T\cos i_T{\sin }^2{u}_T{\cos u}_T}{R_T^6{h}_T^2}\\ 0\\ -\frac{{2h}_T{V}_r}{R_T^3}-\frac{k_J{\sin }^2{i}_T\sin {2u}_T}{R_T^5}\end{array}\right]---\left(6\right)$

式中，R_Lz＝(R_T+x)sin i_Tsin u_T+y sin i_T cos u_T+z cos i_T为洛伦兹航天器地心距在地惯系中Z_I方向的投影，i_T与u_T分别为目标器的轨道倾角与纬度幅角，系数其中地球非球形摄动二阶带谐系数为J₂＝1.0826×10^-3，μ＝3.986005×10¹⁴m³/s²为地球引力常数，R_E为地球半径；

a_C＝[a_R a_S a_W]^T与a_L＝λl＝λ[l_x l_y l_z]^T分别为作用于洛伦兹航天器的控制加速度与洛伦兹加速度，其中λ＝q/m为洛伦兹航天器的荷质比，且

$\left\{\begin{array}{c}{l}_x=\frac{B_0}{R_L^5}{V}_y[({3z}^2-R_L^2)n_z+3z[(R_T+x)n_x+{\mathrm{yn}}_y\left]\right\}-\frac{B_0}{R_L^5}{V}_{z,J_2}\{({3y}^2-R_L^2)n_y+3y[(R_T+x)n_x+{\mathrm{zn}}_z\left]\right\}\\ l_y=\frac{B_0}{R_L^5}{V}_z\left\{\right[3{(R_T+x)}^2-R_L^2]n_z+3(R_T+x)({\mathrm{yn}}_y+{\mathrm{zn}}_z)\}-\frac{B_0}{R_L^5}{V}_x[({3z}^2-R_L^2)n_z+3z[(R_T+x)n_x+{\mathrm{zn}}_z\left]\right\}\\ l_z=\frac{B_0}{R_L^5}{V}_x\{({3y}^2-R_L^2)n_y+3y[(R_T+x)n_x+{\mathrm{zn}}_x\left]\right\}-\frac{B_0}{R_L^5}{V}_y\left\{\right[3{(R_T+x)}^2-R_L^2]n_z+3(R_T+x)({\mathrm{yn}}_y+{\mathrm{zn}}_z)\}\end{array}\right.---\left(7\right)$

式中，B₀＝8.0×10¹⁵T·m³为地磁场地磁矩值。V_rel＝[V_x V_y V_z]^T为洛伦兹航天器与当地磁场相对速度在相对运动坐标系中的描述：

$V_{\mathrm{rel}}=\left[\begin{array}{c}{V}_x\\ V_y\\ V_z\end{array}\right]=\left[\begin{array}{c}{\stackrel{\·}{R}}_T+\stackrel{\·}{x}-y({\mathrm{\ω}}_z-{\mathrm{\ω}}_E\cos i_T)-{\mathrm{z\ω}}_E\sin i_T\cos u_T\\ \stackrel{\·}{y}+(R_T+x)({\mathrm{\ω}}_z-{\mathrm{\ω}}_E\cos i_T)-z({\mathrm{\ω}}_x-{\mathrm{\ω}}_E\sin i_T\sin u_T)\\ \stackrel{\·}{z}+(R_T+x){\mathrm{\ω}}_E\sin i_T\cos u_T+y({\mathrm{\ω}}_x-{\mathrm{\ω}}_E\sin i_T\sin u_T)\end{array}\right]---\left(8\right)$

$\hat{n}={\left[\begin{array}{ccc}{n}_x& n_y& n_z\end{array}\right]}^T$ 为单位磁偶极子矢量在相对运动坐标系中的描述：

$\hat{n}=\left[\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right]=\left[\begin{array}{c}-(\cos \mathrm{\β}\cos u_T+\sin \mathrm{\β}\cos i_T\sin u_T)\sin \mathrm{\α}-{\sin i}_T{\sin u}_T\cos \mathrm{\α}\\ (\cos \mathrm{\β}\sin u_T-\sin \mathrm{\β}\cos i_T\cos u_T)\sin \mathrm{\α}-{\sin i}_T{\cos u}_T\cos \mathrm{\α}\\ \sin \mathrm{\β}\sin i_T\sin \mathrm{\α}-{\cos i}_T\cos \mathrm{\α}\end{array}\right]---\left(9\right)$

式中，β＝Ω_M-Ω_T，其中Ω_M＝ω_Et+Ω₀为地磁轴在地惯系赤道平面内的相位角，ω_E为地球自转角速度，Ω₀为初始时刻地磁轴在地惯系赤道平面内的相位角，Ω_T为目标器轨道升交点赤经，倾斜地磁轴与地球自转轴的夹角为α；

六个轨道参数(R_T,V_r,h_T,i_T,u_T,Ω_T)分别为目标器的轨道半径、径向速度、轨道角动量、轨道倾角、近地点幅角以及升交点赤经，其满足动力学方程

$\left\{\begin{array}{c}{\stackrel{\·}{R}}_T=V_r\\ {\stackrel{\·}{V}}_r=-\mathrm{\μ}/R_T^2+h_T^2/R_T^3-k_J(1-3{\sin }^2{i}_T{\sin }^2{u}_T)/R_T^4\\ {\stackrel{\·}{h}}_T=-\left(k_J{\sin }^2{i}_T\sin 2u_T\right)/R_T^3\\ {\stackrel{\·}{\stackrel{\·}{i}}}_T=-\left(k_J\sin {2i}_T\sin 2u_T\right)/\left({2h}_T{R}_T^3\right)\\ {\stackrel{\·}{u}}_T=h_T/R_T^2+\left({2k}_J\mathrm{co}{s}^2{i}_T\mathrm{si}{n}^2{u}_T\right)/\left(h_T{R}_T^3\right)\\ {\stackrel{\·}{\mathrm{\Ω}}}_T=-\left({2k}_J\cos i_T{\sin }^2{u}_T\right)/\left(h_T{R}_T^3\right)\end{array}\right.---\left(10\right)$

式(2)至式(10)即J₂摄动作用下的洛伦兹航天器非线性相对运动方程；

(2)计算悬停所需的开环控制总加速度：

根据悬停定义，在相对运动坐标系中，追踪器与目标器的相对位置矢量保持不变，因而相对位置各轴分量对时间的各阶导数为零，即

$\left\{\begin{array}{c}\stackrel{\·}{x}=\stackrel{\·}{y}=\stackrel{\·}{z}=0\\ \stackrel{\·\·}{x}=\stackrel{\·\·}{y}=\stackrel{\·\·}{z}=0\end{array}\right.---\left(11\right)$

将式(11)代入式(2)中得，

a_C＝h-λl     (12)

式中，h为实现悬停所需的总控制加速度，即

$h=\left[\begin{array}{c}x({\mathrm{\η}}_L^2-{\mathrm{\ω}}_z^2)-y{\mathrm{\ϵ}}_z+{\mathrm{z\ω}}_x{\mathrm{\ω}}_z+({\mathrm{\ξ}}_L-\mathrm{\ξ})\sin i_T\sin u_T+R_T({\mathrm{\η}}_L^2-{\mathrm{\η}}^2)\\ {\mathrm{x\ϵ}}_z+y({\mathrm{\η}}_L^2-{\mathrm{\ω}}_z^2-{\mathrm{\ω}}_x^2)-z{\mathrm{\ϵ}}_x+({\mathrm{\ξ}}_L-\mathrm{\ξ})\sin i_T\cos u_T\\ {\mathrm{x\ω}}_x{\mathrm{\ω}}_z+{\mathrm{y\ϵ}}_x+z({\mathrm{\η}}_L^2-{\mathrm{\ω}}_x^2)+({\mathrm{\ξ}}_L-\mathrm{\ξ})\cos i_T\end{array}\right]---\left(13\right)$

在步骤四中所述的实现悬停所需的最优洛伦兹航天器荷质比，其计算方法如下：

由于洛伦兹加速度方向与当地磁场方向及洛伦兹航天器相对磁场速度方向垂直，当洛伦兹航天器轨道位置确定，洛伦兹加速度方向也随之确定；若洛伦兹加速度方向与悬停所需加速度方向重合，即l∥h，则若洛伦兹航天器荷质比大小为λ＝||h||/||l||，产生的洛伦兹加速度即可完全抵偿悬停所需加速度；然而，由于洛伦兹力作用方向的局限性，并不能保证洛伦兹加速度方向与悬停所需加速度方向重合，因此，为实现悬停，需施加额外的控制加速度，本发明中假设该额外控制加速度由推力器提供，为最大限度地节省燃料，选取拉格朗日形式的燃料最优目标函数，即：

$J={\∫}_0^{t_f}L[t,\mathrm{\λ}\left(t\right)]\mathrm{dt}={\∫}_0^{t_f}\left|\right|a_C\left|\right|\mathrm{dt}---\left(14\right)$

式中，t_f为悬停任务历时；

求解欧拉-拉格朗日方程：

$\frac{d}{\mathrm{dt}}\left(\frac{\∂L}{\∂\stackrel{\·}{\mathrm{\λ}}}\right)-\frac{\∂L}{\∂\mathrm{\λ}}=0---\left(15\right)$

得到实现燃料最省的洛伦兹航天器荷质比最优轨迹

${\mathrm{\λ}}^{*}\left(t\right)=\left\{\begin{array}{cc}\frac{h\·l}{{\left|\right|l\left|\right|}^2},& \left|\right|l\left|\right|\≠0\\ 0,& \left|\right|l\left|\right|=0\end{array}\right.---\left(16\right)$

当磁场B方向与相对速度V_rel方向平行时，无洛伦兹加速度产生，即||l||＝0，因而，悬停所需加速度将完全由推力器提供，该条件下设定洛伦兹航天器荷质比为0，如式(16)所示；

在步骤五中所述的实现悬停所需最优洛伦兹加速度及推力器提供的最优控制加速度，其计算方法如下：

将洛伦兹航天器最优荷质比λ^*代入洛伦兹加速度定义式，得到实现悬停所需的最优洛伦兹加速度为

$a_L^{*}\left(t\right)=\left\{\begin{array}{cc}\frac{h\·l}{{\left|\right|l\left|\right|}^2}l,& \left|\right|l\left|\right|\≠0\\ 0,& \left|\right|l\left|\right|=0\end{array}\right.---\left(17\right)$

将式(17)代入式(12)，得到最优控制加速度

$a_C^{*}\left(t\right)=\left\{\begin{array}{cc}h-\frac{h\·l}{{\left|\right|l\left|\right|}^2}l,& \left|\right|l\left|\right|\≠0\\ h,& \left|\right|l\left|\right|=0\end{array}\right.---\left(18\right)$

因此，对于带电洛伦兹航天器，实现悬停所需的速度增量为

$\mathrm{\ΔV}={\∫}_0^{t_f}\left|\right|a_C^{*}\left|\right|\mathrm{dt}---\left(19\right)$

而对于非带电航天器，实现悬停所需的控制加速度均由推力器提供，因而所需速度增量为

$\mathrm{\ΔV}={\∫}_0^{t_f}\left|\right|h\left|\right|\mathrm{dt}---\left(20\right)$