近年来，航天器逐渐朝着轻薄化、小型化、模块化发展。微纳卫星(1~50 kg)由于研制周期短、造价低、可批量标准化生产，受到广泛关注。多个微纳卫星通过编队协同工作，实现近地空间以及行星际空间的复杂任务^[1-2]。同时，卫星在轨运行期间出现故障时，通过小卫星对其进行近距离操作，以完成在轨维修、燃料加注、状态监测等工作。对于太空垃圾清除、非合作航天器捕获等任务，也需要近距离操作。为了实现复杂多样的近距离操作任务，交会对接成为关键技术。交会对接一般可以分为4个阶段：远程导引段(15~100 km、近程导引段(0.5~1 km)、最终逼近段(100 m)和对接停靠段。针对交会对接问题，研究成果可以分为2类:一类以转移轨迹设计与优化为目标，使得燃料消耗最小，同时满足推力幅值约束、时间约束等条件^[3-4]; 另一类为相对姿态与位置控制，旨在以较高精度跟踪转移轨迹，同时避免与目标航天器发生碰撞^[5-7]。

目前大部分关于交会对接的研究都针对的是全驱动或者过驱动系统，追踪航天器上安装的推力器可以产生任意方向的控制力，同时姿态控制依靠飞轮和推力器提供力矩。位置控制虽然受姿态的影响，但是在任意姿态下，通过控制分配都可以产生期望推力，即使产生的推力小于期望值，但是可以保证两者方向相同。姿态控制一般不受位置控制影响。

然而对于一些新型卫星，推力方向及大小会受姿态影响。例如太阳帆航天器，其推力需要利用太阳光压产生。不同姿态下，太阳光压不同，从而导致推力不同^[8]。除此之外，用于深空探测、天基发射等任务的微纳卫星，为了降低研制成本、减轻重量、缩小体积，通常只安装一个推力器用于位置控制。这种情况下，位置控制严重依赖于姿态指向，因为姿态指向决定了推力方向。同时，由于推力器只能提供一个方向的推力，在不改变姿态的情况下无法减速，这给交会对接带来了很大挑战。

单推力航天器控制在卫星编队中研究较多。文献[9]针对欠驱动卫星，首先基于势函数设计标称控制力，实现期望构型的同时避免碰撞。然后将标称控制力方向作为期望指向，设计姿态控制器，使得推力器方向与标称控制力重合。这种内-外环控制策略在文献[10]中应用在无人机控制问题上，实现对期望轨迹的跟踪。无人机控制中推力方向与重力相反，因此可以相互抵消，在推力方向的速度可以控制。但是对于卫星，推力方向没有外力与之抵消，因此减速难以实现。Haghighi和Pang^[11]针对微纳卫星编队系统，通过多级控制策略，使得单推力器控制下的微纳卫星形成期望构型。但是这些文献对于末端姿态没有限制，只能实现构型控制。针对悬停任务，文献[12]分别考虑径向及切向无推力的情况，分析了悬停点的范围，并通过反步法实现欠驱动控制。Muralidharan和Emami^[13]设计位置和姿态同步控制器，实现单推力航天器对非合作目标的逼近及绕飞，同时利用这种控制策略完成高轨道向低轨道的转移，可以近似达到燃料最优。

对于绕飞等任务，其最终期望力的方向是朝向目标的，因此可以保证追踪星最终指向目标，整个设计过程只保证位置控制即可。但是对于交会对接，末端的姿态约束取决于对接装置安装位置。要实现以期望的姿态到达指定位置，对于单推力航天器来说具有很大挑战。传统上，对于全推力航天器，完成六自由度控制可以不依赖于轨迹规划，通过直接设计闭环控制律实现位置和姿态镇定。而对于单推力航天器，位置和姿态同时镇定是不可实现的，因为位置的镇定需要控制器在平衡点附近能够提供各个方向的推力，这要求航天器姿态必须变化以提供相应推力。因此，实现六自由度控制必须依靠轨迹规划，且终端速度不能为0。

文献[14]针对单推力航天器近距离操作任务，考虑视线角约束、避碰约束和羽流规避约束，采用脉冲喷气的形式，提出一种改进的快速搜索随机树(RRT^*)轨迹规划策略，以保证轨迹加速度方向与目标方向夹角在一定范围内，从而使得目标可以出现在追踪星的视场中；然后设计有限时间滑模控制器完成姿态重定向，推力方向与期望加速度方向一致。Yoshimura^[15]基于傅里叶级数方法设计了追踪星的运动轨迹，使其以椭圆螺旋线的形式到达绕飞轨迹上，并保证能量最优。同时，保证期望推力方向在惯性系中保持在一定范围内。但是其只考虑了XY平面内的运动情况，并且虽然保证了推力方向，但是是在惯性系中实现的。并不能直接应用在交会对接中。文献[16]针对二维平面内运动的欠驱动系统，依据非完整约束特性设计参考轨迹，并采用滑模控制完成轨迹跟踪。Biggs和Henninger^[17]基于微分几何理论，设计了欠驱动系统的能量最优转移轨迹，给出了轨迹的解析解，通过边界条件，显示求解轨迹参数。该方法可以应用于二维平面交会对接问题。进一步，Henninger和Biggs^[18]通过覆盖映射方法，将六自由度运动方程转化为两组相互耦合的三维方程，然后设计参考转移轨迹，应用于三维空间的交会对接问题。该方法既可以应用于全驱动系统，也可以应用于欠驱动航天器，航天器只在初始时刻进行一次脉冲机动，然后进入推力无控阶段，仅依靠姿态调节完成交会对接。Mitani和Yamakawa^[19]基于(Control Lyapunov Function, CLF)方法，提出一种满足控制力方向约束的制导策略，实现以零速度逼近目标，并可以证明系统稳定性。通过引入满意度函数，使得控制力方向与目标之间夹角小于一个常值。但是追踪星的转移轨迹是以目标为中心的螺旋线。对于交会对接，这种逼近轨迹将导致追踪星与目标发生碰撞。随后，Mitani和Yamakawa在文献[20]中进一步考虑了推力幅值约束。

综上所述，目前对于单推力航天器的交会对接问题研究还有待深入，尤其是针对末端指向固定的情况，难度较大。必须合理设计转移轨迹来实现交会对接。本文针对该问题，首先构建了欠驱动航天器相对位置和姿态动力学方程，然后在不考虑初始推力器方向的情况下采用传统螺旋线设计转移轨迹，满足末端位置和速度要求。进一步考虑初始推力器方向，提出旋转螺旋线设计方法，使其满足初末位置和速度方向要求。最后，设计位置跟踪及姿态跟踪控制器，使得追踪星跟踪上转移轨迹，完成交会对接。

1 航天器相对位置及姿态模型 1.1 相对位置模型

对于交会对接任务，追踪星一般经历调相、近距离交会以及最终接近的过程。为了描述其相对位置，需要建立追踪星与目标星之间的相对位置方程。本文主要研究近距离交会和最终对接阶段，在LVLH坐标系下建立相对位置动力学。LVLH坐标系OX_LY_LZ_L以目标星质心为中心，X_L沿地心指向目标星方向，Z_L沿目标星轨道角速度方向，X_L、Y_L和Z_L成右手法则。交会对接任务及LVLH坐标系示意图如图 1所示。其中r_oc、r_ot、r_tc分别表示地心指向追踪星C、目标星T的矢量以及目标星指向追踪星的矢量；X_cb、X_tb分别为追踪星和目标星的本体系X轴；f为沿追踪星本体X_cb方向的推力矢量。

考虑到整个交会对接过程中，两星距离远小于目标星轨道半径，因此用Clohessy-Wiltshire (C-W)方程不会产生较大误差。两星相对位置方程为^[21]

$ \left\{ {\begin{array}{*{20}{l}} {\ddot x - 2{\omega _0}\dot y - 3\omega _0^2x = \frac{{{f_{{\rm{d}}x}}}}{{{m_{\rm{c}}}}}}\\ {\ddot y + 2{\omega _0}\dot x = \frac{{{f_{{\rm{d}}y}}}}{{{m_{\rm{c}}}}}}\\ {\ddot z + \omega _0^2z = \frac{{{f_{{\rm{d}}z}}}}{{{m_{\rm{c}}}}}} \end{array}} \right. $	（1）

式中：r_tc=[x, y, z]^T为追踪星相对目标星位置；ω₀和r₀为目标星轨道角速度及轨道半径；m_c为追踪星质量；f_d=[f_dx, f_dy, f_dz]^T为LVLH坐标系中追踪星控制力。相对位置动力学可以写成矢量形式：

$ {\mathit{\boldsymbol{\ddot r}}_{{\rm{tc}}}} = \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{r}}_{{\rm{tc}}}} + \mathit{\boldsymbol{N}}{\mathit{\boldsymbol{\dot r}}_{{\rm{tc}}}} + \frac{{{\mathit{\boldsymbol{f}}_{\rm{d}}}}}{{{m_{\rm{c}}}}} $	（2）

式中：

$ \mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} {3\omega _0^2}&0&0\\ 0&0&0\\ 0&0&{ - \omega _0^2} \end{array}} \right];\mathit{\boldsymbol{N}} = \left[ {\begin{array}{*{20}{c}} 0&{2{\omega _0}}&0\\ { - 2{\omega _0}}&0&0\\ 0&0&0 \end{array}} \right] $

1.2 姿态误差描述

与传统交会对接不同，本文中追踪星只安装一个推力器，推力方向沿本体X_cb轴方向。在进行相对位置控制时，为了提供相应控制力，追踪星需要进行姿态调整，将X_cb轴与期望控制力f_d的方向重合。因此，为了方便后续姿态指向控制器设计，需要建立姿态误差方程。定义姿态指向误差为

$ \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{e}}_\theta } = {{[{e_{\theta x}},{e_{{\theta _y}}},{e_{\theta z}}]}^{\rm{T}}}}\\ {{e_{\theta x}} = {\rm{cos}}{\kern 1pt} {\kern 1pt} {\theta _x} - 1 = t_x^{\rm{b}} - 1 = {\mathit{\boldsymbol{R}}_{{\rm{bI}}}}(1,z){\mathit{\boldsymbol{t}}^{\rm{I}}} - 1} \end{array}\\ \begin{array}{*{20}{l}} {{e_{\theta y}} = {\rm{cos}}{\kern 1pt} {\kern 1pt} {\theta _y} = t_y^{\rm{b}} = {\mathit{\boldsymbol{R}}_{{\rm{bI}}}}(2,:){\mathit{\boldsymbol{t}}^{\rm{I}}}}\\ {{e_{\theta z}} = {\rm{cos}}{\kern 1pt} {\kern 1pt} {\theta _z} = t_z^{\rm{b}} = {\mathit{\boldsymbol{R}}_{{\rm{bI}}}}(3,:){\mathit{\boldsymbol{t}}^{\rm{I}}}} \end{array} \end{array} \right. $

（3）

式中：$ \boldsymbol{t} = \frac{{{\boldsymbol{f}_{\rm{d}}}}}{{\parallel {\boldsymbol{f}_{\rm{d}}}\parallel }} $表示期望推力方向；$ {\boldsymbol{t}^{\rm{b}}} = {\left[ {t_x^{\rm{b}}, t_y^{\rm{b}}, t_z^{\rm{b}}} \right]^{\rm{T}}} $和$ {\boldsymbol{t}^{\rm{I}}} = {\left[ {t_x^{\rm{I}}, t_y^{\rm{I}}, t_z^{\rm{I}}} \right]^{\rm{T}}} $分别表示t在追踪星本体系及惯性系的分量；θ_x、θ_y和θ_z分别表示期望推力方向与追踪星本体三轴的夹角，其期望值为θ_xd=0, θ_yd=π/2, θ_zd=π/2;R_bI为追踪星本体系相对惯性系的姿态旋转矩阵，R_bI(j, :)表示R_bI的第j行。根据式(3)，姿态指向误差可以写成

$ {\mathit{\boldsymbol{e}}_\theta } = {\mathit{\boldsymbol{R}}_{{\rm{bI}}}}{\mathit{\boldsymbol{t}}^{\rm{I}}} - {[1,0,0]^{\rm{T}}} $	（4）

对e_θ求导，可得

$ \begin{array}{l} {{\mathit{\boldsymbol{\dot e}}}_\theta } = {{\mathit{\boldsymbol{\dot R}}}_{{\rm{bI}}}}{\mathit{\boldsymbol{t}}^{\rm{I}}} + {\mathit{\boldsymbol{R}}_{{\rm{bI}}}}{{\mathit{\boldsymbol{\dot t}}}^{\rm{I}}} = - \omega _{{\rm{bI}}}^ \times {R_{{\rm{bI}}}}{\mathit{\boldsymbol{t}}^{\rm{I}}} + {R_{{\rm{bI}}}}{{\mathit{\boldsymbol{\dot t}}}^{\rm{I}}} = {\mathit{\boldsymbol{t}}^{\rm{b}}} \times {\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{R}}_{{\rm{bI}}}}{{\mathit{\boldsymbol{\dot t}}}^{\rm{I}}} = \left[ {\begin{array}{*{20}{c}} 0&{ - {e_{\theta z}}}&{{e_{\theta y}}}\\ {{e_{\theta z}}}&0&{ - 1 - {e_{\theta x}}}\\ { - {e_{\theta y}}}&{1 + {e_{\theta x}}}&0 \end{array}} \right]{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{R}}_{{\rm{bI}}}}{{\mathit{\boldsymbol{\dot t}}}^{\rm{I}}} = \mathit{\boldsymbol{F}}({\mathit{\boldsymbol{e}}_\theta }){\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} + {\mathit{\boldsymbol{R}}_{{\rm{bI}}}}{{\mathit{\boldsymbol{\dot t}}}^{\rm{I}}} \end{array} $

（5）

式中：$ \boldsymbol{F}\left( {{\boldsymbol{e}_\theta }} \right) = \left[ {\begin{array}{*{20}{c}} 0&{ - {e_{\theta z}}}&{{e_{\theta y}}}\\ {{e_{\theta z}}}&0&{ - 1 - {e_{\theta x}}}\\ { - {e_{\theta y}}}&{1 + {e_{\theta x}}}&0 \end{array}} \right];{\boldsymbol{\mathop \omega \limits^ \cdot } _{_{{\rm{bI}}}}} $为追踪星本体系相对惯性系的角速度。追踪星转动动力学方程为

$ {\mathit{\boldsymbol{\dot \omega }}_{{\rm{bI}}}} = - \mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^ \times {\mathit{\boldsymbol{J}}_{\rm{s}}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} + \mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}({\mathit{\boldsymbol{\tau }}_{\rm{c}}} + \mathit{\boldsymbol{d}}) $	（6）

式中：J_s为转动惯量；τ_c为控制力矩；d表示外部干扰。

2 交会对接轨迹规划

交会对接任务需要保证追踪星与目标的相对位置和相对姿态。对于三轴都安装有推力器的卫星，其交会对接相对简单，因为姿态与位置控制可以认为是解耦的，通过设计闭环反馈控制，可以实现姿态-位置协同控制，在任意姿态下期望控制力可以由多个推力器产生。但是对于单推力器系统，期望控制力必须依靠姿态转动实现。这种情况下，直接设计闭环反馈控制实现交会对接是极其困难的。因此，需要首先设计位置转移轨迹，再完成轨迹跟踪。

需要说明的是，单推力器系统没有减速机制，即追踪星沿推力器方向的速度不可减小。最终对接时刻要求追踪星到达目标星位置，且本体X_cb轴与LVLH坐标系的Y_L轴重合。因此，最终时刻必须保证相对位置转移轨迹的切线方向与Y_L一致。但是对于位置转移曲线，需要一个向心力完成转弯，除直线以外的任何曲线加速度方向并不能沿着切线方向。这就导致虽然追踪星最终速度方向与Y_L轴重合，但是由于期望加速度方向与Y_L并不重合，本体X_cb轴必须指向期望加速度方向以提供相应的推力。综上所述，在最终对接阶段，相对位置曲线需要设计为直线形式。本节将设计两阶段轨迹转移曲线，第1阶段(近距离交会阶段)以螺旋线形式进行逼近，第2阶段(最终对接阶段)以直线形式进行对接。

2.1 不考虑初始速度方向

首先，不考虑追踪星在LVLH系的初始速度方向，基于螺旋线方程，设计转移曲线为

$ \left\{ {\begin{array}{*{20}{l}} {{x_{\rm{d}}}(t) = r{\rm{cos}}(At + B) + {D_1}}\\ {{y_{\rm{d}}}(t) = r{\rm{sin}}(At + B) + {D_2}}\\ {{z_{\rm{d}}}(t) = {C_0} + {C_1}t} \end{array}} \right. $	（7）

式中：r_tcd(t)=[x_d(t), y_d(t), z_d(t)]^T为LVLH坐标系中追踪星的期望轨迹；t为时间变量，B为参数；C₀和C₁分别为转移轨迹(螺旋线)z方向的初始速度和加速度。式(7)设计的螺旋线在OX_LY_L平面内是一个以r为半径、以(D₁, D₂)为圆心的圆，在Z_L方向以匀速直线运动逼近目标，如图 2所示，其中r_tc^⊥表示与r_tc的垂直矢量；A为从初始点到轨迹切换点之间的角度。

根据最优几何规划理论，式(7)设计的螺旋线可以保证对于给定初始速度方向和末端速度方向的任务，曲线的曲率积分是最小的^[22-24]。对于单推力航天器来说，由于推力方向需要指向曲线加速度方向，如果曲率较大，则航天器需要在很短的时间内完成较大姿态机动，这对姿态控制系统的控制算法和执行机构要求较高，微纳卫星由于成本和体积的限制，很难实现敏捷姿态机动。采用螺旋线轨迹，由于曲率积分最优，可以减小姿态旋转角度。同时，圆形曲线可以通过连续平稳转弯降低期望控制力(方向)的变化，减轻姿态控制负担。为了使追踪星在对接段沿着Y_L方向运动，Z_L方向的速度必须减为0。这就需要在交会段结束之前通过改变Z_L方向的轨迹，使其减速运动。因此，整个交会对接过程分为近距离交会-减速-对接阶段，其轨迹方程为

交会段(0≤t≤t₁)

$ \left\{ {\begin{array}{*{20}{l}} {{x_{\rm{d}}}(t) = r{\rm{cos}}(At + B) + {D_1}}\\ {{y_{\rm{d}}}(t) = r{\rm{sin}}(At + B) + {D_2}}\\ {{z_{\rm{d}}}(t) = {C_0} + {C_1}t} \end{array}} \right. $	（8a）

减速段(t₁＜t≤t₂)

$ \left\{ {\begin{array}{*{20}{l}} {{x_{\rm{d}}}(t) = r{\rm{cos}}(At + B) + {D_1}}\\ {{y_{\rm{d}}}(t) = r{\rm{sin}}(At + B) + {D_2}}\\ {{z_{\rm{d}}}(t) = {C_0} + {C_1}{t_1} + 0.5{a_z}{{(t - {t_1})}^2}} \end{array}} \right. $	（8b）

对接段(t₂＜t≤t_f)

$ \left\{ \begin{array}{l} {x_{\rm{d}}}(t) = 0\\ {y_{\rm{d}}}(t) = r{\rm{sin}}(A{t_2} + B) + {D_2} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} rA{\rm{cos}}(A{t_2} + B)(t - {t_2})\\ {z_{\rm{d}}}(t) = 0 \end{array} \right. $

（8c）

式中:a_z为追踪星延z方向加速度。

以上转移轨迹可以分为2个阶段：第1阶段为交会阶段，追踪星以螺旋线形式对目标逼近。需要说明的是，该阶段包括了减速段，如式(8b)所示。虽然减速段轨迹不是严格的螺旋线，但是其作用时间相对较短，对整体性能影响不大。第2阶段为对接段，追踪星以匀速直线运动逼近目标，完成最终对接。

给定初始位置(x₀, y₀, z₀)、末端位置(x_f, y_f, z_f)、末端速度方向v_f、对接起始位置(轨迹切换点)y₂₀以及各阶段起始结束时刻t₁、t₂，可以解算出轨迹参数r、A、B、D₁、D₂、C₀、C₁和a_z为

$ \left\{ {\begin{array}{*{20}{l}} {{D_2} = {y_{20}},{D_1} = \frac{{x_0^2 + y_0^2 - 2{y_0}{y_{20}} + y_{20}^2}}{{2{x_0}}}}\\ {r = \sqrt {{{({x_0} - {D_1})}^2} + {{({y_0} - {D_2})}^2}} }\\ {B = {\rm{arctan}}\frac{{{y_0} - {D_2}}}{{{x_0} - {D_1}}}}\\ {A = \frac{{{k_2}\pi - B}}{{{t_2}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {k_2} = 0, \pm 1, \pm 2, \cdots }\\ {{C_0} = {z_0},{C_1} = \frac{{ - 2{C_0}}}{{{t_1} + {t_2}}}}\\ {{a_z} = \frac{{ - {C_1}}}{{{t_2} - {t_1}}}} \end{array}} \right. $

（9）

需要注意的是，k₂的目的是保证D₁A＜0。否则最终的速度方向为-Y_L。根据以上求得的参数，可以得到终端时刻

$ {t_f} = {t_2} + \frac{{{y_{20}}}}{{rA{\rm{cos}}(A{t_2} + B)}} $	（10）

虽然这种情况下曲线形状唯一，但是由于终端时刻没有限制，因此曲线的速度是可以通过调节A参数改变的，A决定了圆周运动的角速度，其并不影响曲线形状，在曲线形状不变的情况下，A越大，t_f越小。若限定终端时刻t_f，则

$ A = \pm \frac{{|{y_{20}}|}}{{r({t_{\rm{f}}} - {t_2})}} $	（11）

其正负号选择根据条件D₁A＜0。

另一种情况，如果限定t_f、t₁和t₂，而对接起始位置y₂₀自由，则同样也可以得到一组解集:

$ \left\{ {\begin{array}{*{20}{l}} {{{({x_0} - {D_1})}^2} + {{({y_0} - {D_2})}^2} = D_1^2}\\ {r = |{D_1}|}\\ {B = {\rm{arctan}}\frac{{{y_0} - {D_2}}}{{{x_0} - {D_1}}},A = \pm \frac{{{D_2}}}{{{D_1}({t_{\rm{f}}} - {t_2})}}}\\ {{C_0} = {z_0},{C_1} = \frac{{2{C_0}}}{{{t_2} - 3{t_1}}}}\\ {{a_z} = \frac{{ - {C_1}}}{{{t_2} - {t_1}}}} \end{array}} \right. $	（12）

此时曲线形状不唯一，可以给定任意圆心位置和半径。由于|D₂|决定了最后对接段的速度，因此其值越小，最终追踪星与目标星的相对速度越小。

注1   螺旋线转移轨迹虽然可以保证曲率积分最优，但是由于没有考虑相对位置动力学，其难以保证燃料消耗的最优性。与传统的交会轨迹优化方法相比，螺旋线的设计方法更加侧重曲线的几何特性而非动力学系统。因此，基于螺旋线设计的转移轨迹虽然有利于姿态控制，但是对于相对位置控制优势不明显。

2.2 考虑初始速度方向

2.1节中，虽然给出了螺旋线形式的转移轨迹，但是并没有考虑追踪星初始X_cb方向。为了产生转移轨迹切向的控制力，追踪星在初始时刻需要进行姿态机动。此时，转移轨迹曲率积分最小的优越性将失去意义，因为追踪星需要跟踪上转移轨迹，而对于只安装有小推力、单方向推力器的追踪星，完成轨迹跟踪有较大挑战性，其需要较长时间从初始位置和速度转移到转移轨迹上。

因此，考虑追踪星推力器初始指向是十分必要的。本节仍然采用螺旋线+直线形式的转移轨迹，其中直线轨迹不施加控制，在轨迹切换后以固定速度逼近目标，如图 3所示。

本节在交会段所考虑的边界条件如下：

初末位置约束：

$ \left\{ {\begin{array}{*{20}{l}} {{x_{\rm{d}}}(0) = {x_0}}\\ {{y_{\rm{d}}}(0) = {y_0},}\\ {{z_{\rm{d}}}(0) = {z_0}} \end{array}\left\{ {\begin{array}{*{20}{l}} {{x_{\rm{d}}}({t_1}) = 0}\\ {{y_{\rm{d}}}({t_1}) = {y_{20}}}\\ {{z_{\rm{d}}}({t_1}) = 0} \end{array}} \right.} \right. $	（13a）

初末速度方向约束：

$ \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} {[{{\dot x}_{\rm{d}}}(0),{{\dot y}_{\rm{d}}}(0),{{\dot z}_{\rm{d}}}(0)]^{\rm{T}}} = {k_{{\rm{v1}}}}{\mathit{\boldsymbol{x}}_{{\rm{cb}}}}\\ {\left[ {{{\dot x}_{\rm{d}}}({t_1}),{{\dot y}_{\rm{d}}}({t_1}),{{\dot z}_{\rm{d}}}({t_1})} \right]^{\rm{T}}} = \end{array}\\ {{k_{{\rm{v2}}}}{\mathit{\boldsymbol{v}}_{\rm{f}}} = {k_{{\rm{v2}}}}{{[0,1,0]}^{\rm{T}}}} \end{array}} \right. $

（13b）

式中：y₂₀为最后对接段开始位置，即轨迹切换点：x_cb为追踪星本体X_b轴的单位矢量。在t₁时刻转移轨迹从螺旋线形式切换为直线。k_v1和k_v2为两个常值系数，其中k_v1表示转移轨迹在追踪星初始位置点的导数(速度)幅值；k_v2表示对接段起始时刻，追踪星相对目标的速度大小，该值越大，对接段持续时间越短。

在式(13)的边界条件下，如果仍然利用2.1节所设计的常规螺旋线设计转移轨迹，无法获得满足所有约束条件的解。以X_LOY_L平面运动为例，当限定轨迹初末时刻的切线方向后，圆心位置及半径就可以唯一确定，然而圆心位置同时还要保证在追踪星与目标连线的垂直平分线上，这两个条件一般情况下是矛盾的，导致曲线参数无解。本节将对传统螺旋线进行姿态旋转，在不改变曲线空间形状的前提下使其满足式(13)的边界条件。

对于交会阶段，设计转移轨迹为

$ \left[ {\begin{array}{*{20}{l}} {{x_{\rm{d}}}(t)}\\ {{y_{\rm{d}}}(t)}\\ {{z_{\rm{d}}}(t)} \end{array}} \right] = \mathit{\boldsymbol{R}}(\psi ,\theta ,\varphi )\left[ {\begin{array}{*{20}{c}} {r{\rm{cos}}(At + B) + {D_1}}\\ {r{\rm{sin}}(At + B) + {D_2}}\\ {{C_0} + {C_1}t + {C_2}{t^2}} \end{array}} \right] $	（14）

式中：C_z为转移轨迹(螺旋线)z方向初始加速度；R(ψ, θ, φ)=R(ψ)R(θ)R(φ)为姿态旋转矩阵，具体表达式为

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{R}}(\varphi ) = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&{{\rm{cos}}\varphi }&{{\rm{sin}}\varphi }\\ 0&{ - {\rm{sin}}\varphi }&{{\rm{cos}}\varphi } \end{array}} \right]}\\ {\mathit{\boldsymbol{R}}(\theta ) = \left[ {\begin{array}{*{20}{c}} {{\rm{cos}}\theta }&0&{ - {\rm{sin}}\theta }\\ 0&1&0\\ {{\rm{sin}}\theta }&0&{{\rm{cos}}\theta } \end{array}} \right]}\\ {\mathit{\boldsymbol{R}}(\psi ) = \left[ {\begin{array}{*{20}{c}} {{\rm{cos}}\psi }&{{\rm{sin}}\psi }&0\\ { - {\rm{sin}}\psi }&{{\rm{cos}}\psi }&0\\ 0&0&1 \end{array}} \right]} \end{array}} \right. $

（15）

式中:φ、θ和ψ分别为绕X、Y和Z轴欧拉转角。

通过将LVLH中以Z_L为中心的螺旋线进行1-2-3转序的欧拉旋转，改变了中心线的位置。由于最后一次旋转是绕Z轴(新的中心线)旋转，不改变螺旋线的位置，可以通过增加z_d(t)方向的自由度达到同样的效果。因此，只采用两次欧拉旋转即可，令ψ=0。同时这样也可以降低计算复杂度。此时，转移轨迹中的未知参数为r、A、B、D₁、D₂、C₀、C₁、C₂、φ、θ以及k_v1和k_v2。由于涉及复杂的三角函数，无法得到相关参数的解析解，只能利用相关数学软件求得数值解。

将空间螺旋线进行旋转后，螺旋线的中心轴在LVLH坐标系中发生了转动，这使得螺旋线在X_LOY_L平面内的投影为不规则曲线，增加了曲线自由度，从而可以满足更复杂的边界条件。需要说明的是，经过旋转后的曲线在空间中依然保持了螺旋线的形状，因此仍然具有曲率积分最小的性质。旋转螺旋线如图 4所示，其中传统螺旋线(实线)以Z_L为中心，经过不同欧拉旋转后，得到一系列旋转螺旋线(虚线)。

根据式(13)和式(14)，可得

$ \left\{ \begin{array}{l} {t_1} + \left| {\frac{{{y_{20}}}}{{{v_y}({t_1})}}} \right| = {t_{\rm{f}}}\\ {v_y}({t_1}) = rA{\kern 1pt} {\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} \varphi {\kern 1pt} {\kern 1pt} {\rm{cos}}(A{t_1} + B) + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({C_1} + 2{C_2}{t_1}){\rm{sin}}{\kern 1pt} {\kern 1pt} \varphi \end{array} \right. $

（16）

式中:v_y(t₁)为转移轨迹在t₁时刻的Y方向速度。

式(16)说明当终端时刻t_f给定后，轨迹切换时刻t₁与最终逼近段距离y₂₀之间的关系。为了给出在交会阶段转移轨迹所需要的控制力，将式(14)代入相对位置动力学中，从而得到

$ \left[ {\begin{array}{*{20}{l}} {f_{{\rm{d}}x}^\prime }\\ {f_{{\rm{d}}y}^\prime }\\ {f_{{\rm{d}}z}^\prime } \end{array}} \right] = \mathit{\boldsymbol{R}}(\psi ,\theta ,\varphi )\left[ {\begin{array}{*{20}{l}} {{f_{{\rm{d}}x}}}\\ {{f_{{\rm{d}}y}}}\\ {{f_{{\rm{d}}z}}} \end{array}} \right] $	（17）

式中：f_dx、f_dy、f_dz表示未经旋转曲线所需控制力，具体表示为

$ \left\{ {\begin{array}{*{20}{l}} {{f_{{\rm{d}}x}} = - {m_{\rm{c}}}(r{A^2} + 2{\omega _0}rA + 3\omega _0^2r){\rm{cos}}(At + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} B) - {m_{\rm{c}}}3\omega _0^2{D_1}}\\ {{f_{{\rm{d}}y}} = - {m_{\rm{c}}}(r{A^2} + 2{\omega _0}rA){\rm{sin}}(At + B)}\\ {{f_{{\rm{d}}z}} = 2{m_{\rm{c}}}{C_2} + {m_{\rm{c}}}\omega _0^2({C_0} + {C_1}t + {C_2}{t^2})} \end{array}} \right. $

（18）

注2  为了使得方程有解，在采用旋转螺旋线时将z_d(t)中的C₀+C₁t变为了C₀+C₁t+C₂t²，虽然这不是严格意义上曲率积分最小的螺旋线，但是由于在近距离交会阶段追踪星和目标星基本在同一轨道平面内，两者距离矢量在Z_L上的分量很小，z_d(t)并不起决定作用，因此在一定程度上仍然能够保证曲线的曲率积分最优性。

3 轨迹跟踪控制

为了能够跟踪上转移轨迹，需要设计相对位置控制器。同时，由于只有在本体X_cb方向可以产生推力，必须通过姿态控制使得本体X_cb轴与期望控制力方向重合。

3.1 位置跟踪控制

定义跟踪误差e_x=r_tc－r_tcd，r_tcd为式(7)和式(14)设计的轨迹。根据式(2)的相对位置动力学，得到位置跟踪误差动力学为

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot e}}}_x} = {\mathit{\boldsymbol{e}}_v}}\\ {{{\mathit{\boldsymbol{\dot e}}}_v} = \mathit{\boldsymbol{M}}({\mathit{\boldsymbol{r}}_{{\rm{tc}}}}){\mathit{\boldsymbol{r}}_{{\rm{tc}}}} + \mathit{\boldsymbol{N}}{{\mathit{\boldsymbol{\dot r}}}_{{\rm{tc}}}} + \frac{{{\mathit{\boldsymbol{f}}_{\rm{d}}}}}{{{m_{\rm{c}}}}} - {{\mathit{\boldsymbol{\ddot r}}}_{{\rm{tcd}}}}} \end{array}} \right. $

（19）

设计PD控制器为

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{f}}_{\rm{d}}} = - {m_{\rm{c}}}(\mathit{\boldsymbol{M}}({\mathit{\boldsymbol{r}}_{{\rm{tc}}}}){\mathit{\boldsymbol{r}}_{{\rm{tc}}}} + \mathit{\boldsymbol{N}}{{\mathit{\boldsymbol{\dot r}}}_{{\rm{tc}}}}) + {m_{\rm{c}}}{{\mathit{\boldsymbol{\ddot r}}}_{{\rm{tcd}}}} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {m_{\rm{c}}}({k_{\rm{p}}}{\mathit{\boldsymbol{e}}_{\rm{x}}} + {k_{\rm{d}}}{\mathit{\boldsymbol{e}}_v})} \end{array} $

（20）

式中:e_v为速度跟踪误差; k_p和k_d分别为位置跟踪控制器比例项和阻尼项系数。

在此控制器下，跟踪误差可渐近收敛为0。但是上述控制器是在LVLH坐标系中设计的，需要将其转化为本体系中的控制力。考虑到姿态控制虽然是快回路，其仍然需要一定时间完成指向控制，追踪星并不能实时产生如式(20)的控制信号。设计本体系中的推力为^[10]

$ {\mathit{\boldsymbol{f}}^{\rm{b}}} = {\left[ {\begin{array}{*{20}{l}} {c({\theta _x})\left\| {{\mathit{\boldsymbol{f}}_{\rm{d}}}} \right\|,0,0} \end{array}} \right]^{\rm{T}}} $	（21）

式中：$ c\left( {{\theta _x}} \right) = \frac{{l - \left( {1 - \cos {\theta _x}} \right)}}{l} $; θ_x表示追踪星本体X_cb轴与期望推力方向f_d的夹角；l为>2的常数；θ_x较大时，表示追踪星无法提供期望推力，此时如果直接在X_cb上施加||f_d||，则会严重影响控制性能。当X_cb与期望推力方向重合后，c(θ_x)=1。通过c(θ_x)，充分考虑了姿态跟踪误差的影响，可以有效提升控制效果。

3.2 姿态跟踪控制

式(20)设计的位置跟踪控制器为姿态控制提供了参考信号，姿态控制的目标就是使得X_cb轴与f_d重合。这本质上是姿态指向控制问题^[25]。本节对传统CLF进行改进^[26-28]，设计基于CLF的滑模控制律，同时保证姿态指向控制的最优性与鲁棒性。

将姿态跟踪误差动力学方程式(5)和式(6)写成如式(22)的仿射形式：

$ \mathit{\boldsymbol{\dot x}} = \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}}) + \mathit{\boldsymbol{g}}(\mathit{\boldsymbol{x}})\mathit{\boldsymbol{u}} $	（22）

式中：

$ \mathit{\boldsymbol{x}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_\theta }}\\ {{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}}} \end{array}} \right],\mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}}) = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}}({\mathit{\boldsymbol{e}}_\theta }){\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} + {\mathit{\boldsymbol{R}}_{{\rm{bI}}}}{{\mathit{\boldsymbol{\dot t}}}^{\rm{I}}}}\\ { - \mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^ \times {\mathit{\boldsymbol{J}}_{\rm{s}}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}}} \end{array}} \right] $

$ \mathit{\boldsymbol{g}} = \left[ {\begin{array}{*{20}{c}} {{{\bf{0}}_{3 \times 3}}}\\ {\mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}} \end{array}} \right],\mathit{\boldsymbol{u}} = {\mathit{\boldsymbol{\tau }}_{\rm{c}}} $

引理1  定义如式(23)的Lyapunov函数

$ \mathit{\boldsymbol{V}} = \frac{a}{2}\mathit{\boldsymbol{e}}_\theta ^{\rm{T}}{\mathit{\boldsymbol{e}}_\theta } + b\mathit{\boldsymbol{e}}_\theta ^{\rm{T}}\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{J}}_{\rm{s}}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} + \frac{c}{2}\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}}{\mathit{\boldsymbol{J}}_{\rm{s}}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} $	（23）

如果a>0，b>0，c>0，则V是一个CLF^[29]。

证明：

由于, $ {\boldsymbol{F}^{\rm{T}}}{\boldsymbol{e}_\theta } = \boldsymbol{P}{\boldsymbol{e}_\theta }, \boldsymbol{P} = \left[ {\begin{array}{*{20}{c}} 0&0&0\\ 0&0&1\\ 0&{ - 1}&0 \end{array}} \right] $, V可以写成

$ \mathit{\boldsymbol{V}} = \left[ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{e}}_\theta ^{\rm{T}}}&{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\frac{a}{2}\mathit{\boldsymbol{I}}}&{\frac{b}{2}{\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{J}}_{\rm{s}}}}\\ {\frac{b}{2}{\mathit{\boldsymbol{J}}_{\rm{s}}}\mathit{\boldsymbol{P}}}&{\frac{c}{2}{\mathit{\boldsymbol{J}}_{\rm{s}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{e}}_\theta }}\\ {{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}}} \end{array}} \right] $

（24）

如果，$ \frac{{ac}}{4}{\boldsymbol{J}_{\rm{s}}} > \frac{{{b^2}}}{4}{\boldsymbol{J}_{\rm{s}}}{\boldsymbol{P}^{\rm{T}}}\boldsymbol{P}{\boldsymbol{J}_{\rm{s}}} $则V>0。下面将证明L_gV=0⇒L_fV＜0。V对x的偏导数为

$ \frac{{\partial {\mathit{\boldsymbol{V}}^{\rm{T}}}}}{{\partial \mathit{\boldsymbol{x}}}} = [a\mathit{\boldsymbol{e}}_\theta ^{\rm{T}} + b\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}}{\mathit{\boldsymbol{J}}_{\rm{s}}}\mathit{\boldsymbol{P}}\quad b\mathit{\boldsymbol{e}}_\theta ^{\rm{T}}\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{J}}_{\rm{s}}} + c\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}}{\mathit{\boldsymbol{J}}_{\rm{s}}}] $

（25）

从而得到

$ \left\{ {\begin{array}{*{20}{l}} {{L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}} = \frac{{\partial {\mathit{\boldsymbol{V}}^{\rm{T}}}}}{{\partial \mathit{\boldsymbol{x}}}}\mathit{\boldsymbol{g}} = b\mathit{\boldsymbol{e}}_\theta ^{\rm{T}}\mathit{\boldsymbol{F}} + c\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}}}\\ {{L_\mathit{\boldsymbol{f}}}\mathit{\boldsymbol{V}} = \frac{{\partial {\mathit{\boldsymbol{V}}^{\rm{T}}}}}{{\partial \mathit{\boldsymbol{x}}}}\mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}}) = (a\mathit{\boldsymbol{e}}_\theta ^{\rm{T}} + b\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}}{\mathit{\boldsymbol{J}}_{\rm{s}}}\mathit{\boldsymbol{P}})\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (b\mathit{\boldsymbol{e}}_\theta ^{\rm{T}}\mathit{\boldsymbol{F}} + c\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}})\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^ \times {\mathit{\boldsymbol{J}}_{\rm{s}}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}}} \end{array}} \right. $

（26）

如果L_gV=0，则有

$ b\mathit{\boldsymbol{e}}_\theta ^{\rm{T}}\mathit{\boldsymbol{F}} + c\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}} = b\mathit{\boldsymbol{e}}_\theta ^{\rm{T}}{\mathit{\boldsymbol{P}}^{\rm{T}}} + c\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}} = {{\bf{0}}_{1 \times 3}} $	（27）

继而

$ \begin{array}{l} {L_\mathit{\boldsymbol{f}}}\mathit{\boldsymbol{V}} = (a\mathit{\boldsymbol{e}}_\theta ^{\rm{T}} + b\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}}{\mathit{\boldsymbol{J}}_{\rm{s}}}\mathit{\boldsymbol{P}})\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} = - \frac{{ac}}{b}\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {b\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^{\rm{T}}{\mathit{\boldsymbol{J}}_{\rm{s}}}\mathit{\boldsymbol{PF}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} \le - \frac{{ac}}{b}{{\left\| {{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}}} \right\|}^2} + b{{\left\| {{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}}} \right\|}^2} \cdot }\\ {\left\| {{\mathit{\boldsymbol{J}}_{\rm{s}}}\mathit{\boldsymbol{PF}}} \right\| = \left( { - \frac{{ac}}{b} + b\sqrt {{\lambda _{{\rm{max}}}}({\mathit{\boldsymbol{J}}_{\rm{s}}}\mathit{\boldsymbol{PF}}{\mathit{\boldsymbol{F}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{J}}_{\rm{s}}})} } \right) \cdot }\\ {{{\left\| {{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}}} \right\|}^2}} \end{array} \end{array} $

（28）

记λ_max(J_SPFF^TP^TJ_S)=η，对于∀ω_bI≠0如果ac>b²η，则L_fV＜0。

根据CLF理论^[30]，得到CLF控制器为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}^*} = \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{u}}_1^*}&{{L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}} \ne {\bf{0}}}\\ {\mathit{\boldsymbol{u}}_2^*}&{{L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}} = {\bf{0}}} \end{array}} \right.}\\ {\mathit{\boldsymbol{u}}_1^* = - {\mathit{\boldsymbol{R}}^{ - 1}} \cdot }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{L_\mathit{\boldsymbol{f}}}\mathit{\boldsymbol{V}} + \sqrt {{{({L_\mathit{\boldsymbol{f}}}\mathit{\boldsymbol{V}})}^2} + l(\mathit{\boldsymbol{x}})({L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}}){\mathit{\boldsymbol{R}}^{ - 1}}{{({L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}})}^{\rm{T}}}} }}{{({L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}}){\mathit{\boldsymbol{R}}^{ - 1}}{{({L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}})}^{\rm{T}}}}} \cdot }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{({L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}})}^{\rm{T}}}}\\ {\mathit{\boldsymbol{u}}_2^* = {\bf{0}}} \end{array}} \right. $

（29）

l(x)=x^TQx，指标函数

$ J = \int_0^\infty {(l(} \mathit{\boldsymbol{x}}) + {\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{Ru}}){\rm{d}}t $	（30）

式中:R和Q为二次型指标函数权值矩阵; u为姿态跟踪控制输入。

传统CLF控制中，当L_gV=0时，为了避免奇异，控制切换为0。而L_gV=0实际上可以理解为一个线性滑模面，即

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{s}} = b\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{e}}_\theta } + c{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} = {\bf{0}}}\\ {{s_1}:{\omega _x} = 0}\\ {{s_2}:b{e_{\theta z}} + c{\omega _y} = 0}\\ {{s_3}: - b{e_{\theta y}} + c{\omega _z} = 0} \end{array}} \right. $	（31）

由于传统CLF控制对于干扰的鲁棒性较差，只能被动抑制干扰，缺乏主动抑制干扰及不确定性的机制。而滑模控制通过不断切换控制信号，将系统误差限制在一定区域内，有效提升系统的控制精度^[31-32]。因此，将滑模控制与CLF结合，设计控制器为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}^*} = \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{u}}_1^*}&{\left\| {{L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}}} \right\| \ge {\delta _1}}\\ {\mathit{\boldsymbol{u}}_2^*}&{\left\| {{L_\mathit{\boldsymbol{g}}}\mathit{\boldsymbol{V}}} \right\| < {\delta _1}} \end{array}} \right.}\\ {\mathit{\boldsymbol{u}}_2^* = \mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^ \times {\mathit{\boldsymbol{J}}_{\rm{s}}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} - \frac{b}{c}{\mathit{\boldsymbol{J}}_{\rm{s}}}\mathit{\boldsymbol{PF}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} - \frac{{{k_1}}}{c}{\mathit{\boldsymbol{J}}_{\rm{s}}} {\rm{sig}} {{(\mathit{\boldsymbol{s}})}^{\frac{m}{n}}}} \end{array}} \right. $

（32）

式中：u₁^*为传统CLF控制，如式(29)中所示；u₂^*中用$ {\rm{sig}}{\left( \boldsymbol{s} \right)^{\frac{m}{n}}} $代替符号函数，目的是使得控制信号变得平滑；m和n为奇数，满足0＜m＜n；δ₁为控制切换条件。

定理1  考虑外部干扰力矩d，在CLF控制律u₁^*的作用下，如果系统状态在t_s时刻进入区域Ω₁，则通过滑模控制律u₂^*，滑模变量将始终保持在Ω₁，并在有限时间T内收敛至Ω₂⊂Ω₁，当时间t→∞时，姿态误差变量e_θ将收敛至域Ω₃。式中：

$ \left\{ \begin{array}{l} {\varOmega _1}:\{ \mathit{\boldsymbol{s}}|\left\| \mathit{\boldsymbol{s}} \right\| \le {\delta _1}\} \\ {\varOmega _2}:\{ \mathit{\boldsymbol{s}}|\left\| \mathit{\boldsymbol{s}} \right\| < {\delta _2}\} ,{\delta _2} = {[\kappa k_1^{ - 1}{(1 - \varepsilon )^{ - 1}}]^{\frac{n}{{m + n}}}},\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \kappa = c\bar d{\delta _1}\left\| {\mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}} \right\|\\ {\varOmega _3}:\left\{ {{\mathit{\boldsymbol{e}}_\theta }|\left\| {\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{e}}_\theta }} \right\| < \frac{{c{\delta _2}}}{b}} \right\} \end{array} \right. $

（33）

证明：

步骤1

构造如下Lyapunov函数并求导

$ {V_1} = \frac{1}{2}{\mathit{\boldsymbol{s}}^{\rm{T}}}\mathit{\boldsymbol{s}} $	（34）

$ \begin{array}{l} \begin{array}{*{20}{l}} {{{\dot V}_1} = {\mathit{\boldsymbol{s}}^{\rm{T}}}[b\mathit{\boldsymbol{PF\omega }} + c\mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}( - \mathit{\boldsymbol{\omega }}_{{\rm{bI}}}^ \times {\mathit{\boldsymbol{J}}_{\rm{s}}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} + \mathit{\boldsymbol{u}}_2^* + \mathit{\boldsymbol{d}})] = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{s}}^{\rm{T}}}( - {k_1} {\rm{sig}} {{(\mathit{\boldsymbol{s}})}^{\frac{m}{n}}} + c\mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}\mathit{\boldsymbol{d}}) = - {k_1}\sum\limits_{i = 1}^3 | {s_i}{|^{\frac{{m + n}}{n}}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} c{\mathit{\boldsymbol{s}}^{\rm{T}}}\mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}\mathit{\boldsymbol{d}} \le - {k_1}{{\left\| \mathit{\boldsymbol{s}} \right\|}^{\frac{{m + n}}{n}}} + c\bar d\left\| \mathit{\boldsymbol{s}} \right\|\left\| {\mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}} \right\| = } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - \left\| \mathit{\boldsymbol{s}} \right\|({k_1}{\left\| \mathit{\boldsymbol{s}} \right\|^{\frac{m}{n}}} - c\bar d\left\| {\mathit{\boldsymbol{J}}_{\rm{s}}^{ - 1}} \right\|) \end{array} $

（35）

式中:d为干扰上界。滑模变量s将收敛至$ \parallel \boldsymbol{s}\parallel ＜ \gamma = {\left( {\frac{{c\parallel \boldsymbol{J}_s^{ - 1}\parallel \overline d }}{{{k_1}}}} \right)^{\frac{n}{m}}}, 记 \parallel \boldsymbol{s}\parallel \le \parallel {\boldsymbol{s}_{t = {t_{\rm{s}}}}}\parallel = {\delta _1} $，则式(35)可以写为

$ {\dot V_1} \le - {k_1}{\left\| \mathit{\boldsymbol{s}} \right\|^{\frac{{m + n}}{n}}} + c\bar d{\delta _1}\left\| {\left. {J_{\rm{s}}^{ - 1}} \right|} \right. = - {k_1}{2^{\frac{{m + n}}{{2n}}}}V_1^{\frac{{m + n}}{{2n}}} + \kappa $	（36）

式中：κ=cdδ₁||J_s^-1||为常数。根据实际有限时间理论，s在有限时间$ T \le \frac{{n\delta _1^{\frac{{n - m}}{n}}}}{{{k_1}\varepsilon \left( {n - m} \right)}} $内收敛至区域Ω₂，ε∈(0, 1)。

步骤2

构造V₂如下：

$ {V_2} = \frac{1}{2}\mathit{\boldsymbol{e}}_\theta ^{\rm{T}}{\mathit{\boldsymbol{e}}_\theta } $	（37）

对其求导得

$ \begin{array}{*{20}{l}} {{{\dot V}_2} = \mathit{\boldsymbol{e}}_\theta ^{\rm{T}}\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{\omega }}_{{\rm{bI}}}} = \mathit{\boldsymbol{e}}_\theta ^{\rm{T}}\mathit{\boldsymbol{P}}\left( {\mathit{\boldsymbol{s}} - \frac{b}{c}\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{e}}_\theta }} \right) \le }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\| {\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{e}}_\theta }} \right\|\left( {{\delta _2} - \frac{b}{c}\left\| {\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{e}}_\theta }} \right\|} \right)} \end{array} $

（38）

因此，$ \mathop {\lim }\limits_{t \to \infty } \parallel \boldsymbol{P}{\boldsymbol{e}_\theta }\parallel \le \frac{{c{\delta _2}}}{b} $，指向误差e_θ将收敛至Ω₃。

4 仿真分析

本节将对本文提出的轨迹规划和跟踪算法进行仿真验证。追踪星参数为^[2]

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{J}}_{\rm{s}}} = 0.01\left[ {\begin{array}{*{20}{c}} {100.9}&0&0\\ 0&{25.1}&0\\ 0&0&{91.6} \end{array}} \right]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{kg}} \cdot {{\rm{m}}^{\rm{2}}}}\\ {{m_{\rm{c}}} = 22.82{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{kg}}} \end{array}} \right. $

追踪星与目标星轨道参数如表 1所示。追踪星轨迹跟踪控制器参数如表 2所示。环境干扰力矩为

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{d}} = {{10}^{ - 4}} \times [{\rm{sin}}(0.001t) - {\rm{cos}}(0.001t)}\\ {{\rm{sin}}(0.001t){]^{\rm{T}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{N}} \cdot {\rm{m}}} \end{array} $

仿真起始时刻令追踪星X_cb轴在LVLH系中的指向为x_cb=[0.1130.991 4－0.066 3]^T。

4.1 不考虑初始速度方向

首先，不考虑转移轨迹的初始方向，利用2.1节中的方法设计参考轨迹。令t₁=8 000 s，t₂=8 200 s，最终逼近段开始位置为(0, -20, 0) m。根据式(10)，可以得到终端时刻t_f=8 311 s。

如图 5所示，实线为式(8)设计的参考轨迹，其分为3段，第1段为近距离交会段(0~8 000 s), 该阶段参考轨迹以螺旋线的形式朝目标运动。为了让Z_L方向的速度减为0，从8 000 s开始进入减速段，直至t=8 200 s，此时参考轨迹到达Y_L轴，且速度方向沿Y_L轴正向。最后，参考轨迹为一直线，匀速逼近目标，消耗时间为111 s。由于在设计参考轨迹时未考虑追踪星初始速度及初始指向，追踪星为了跟踪参考轨迹，需要调整姿态，施加相应控制力进行位置控制。对于微纳卫星，其推力器最大输出幅值一般较小，因此需要较长时间来跟踪上转移轨迹。从图中可以看出，最终追踪星到达目标星位置，且X_cb指向Y_L，符合对接姿态要求。

图 6和图 7分别给出了追踪星对参考轨迹的跟踪误差以及相对目标的位置。可以看出，追踪星需要4 000 s来完成轨迹跟踪。图 7局部放大图是最后逼近段的相对位置变化，可以看出，在初始Y_L方向速度下，追踪星保持该速度驶向目标，最终到达目标附近。需要说明的是，最终产生大约为0.5 m的位置误差，一方面是由于跟踪误差导致的，而跟踪误差很大程度上是由于姿态跟踪误差导致，姿态误差会导致追踪星无法提供期望控制力；令一方面，在最终逼近段中，沿着X_L方向会有外力作用在追踪星上，而为了达到最终对接姿态，追踪星此时推力方向只能沿Y_L, 从而导致X_L方向的误差。如果为了降低最终位置误差，可以令逼近段起始点距离最终对接位置近一些，即减小|y₂₀|。但是这样会减小最终姿态调整的时间，有可能导致追踪星达到对接位置后姿态仍不满足对接要求。图 8为追踪星在LVLH坐标系中的实际速度以及转移轨迹的速度。在逼近段，追踪星距离目标较近，此时如果为了减速而将推力器朝向目标，不能保证对接姿态要求，同时推力器羽流也会污染目标星。因此，追踪星在逼近段不能减速，最终以0.1 m/s的速度与目标星完成对接。图 9是转移轨迹所需要的推力，将转移轨迹代入式(2)即可获得。所需推力幅值很大程度上取决于任务要求时间，即t₁及t_f。本仿真算例中，追踪星需要在8 311 s时间内连续提供约10 mN的推力，这对于现有的化学工质推力器来说是无法实现的，必须依靠电推力器或冷气推进器。电推力器的比冲高、体积小、推力连续可调，可长时间提供小推力，适合用在立方体卫星轨道转移或保持任务中^[33]。实际任务中，可以根据不同推力器特性，调整t₁及t_f，以使得推力器提供的推力能够跟踪上转移轨迹。

为了跟踪参考轨迹，图 10给出了根据式(20)设计的位置控制力曲线，表示在LVLH系中。在初始时刻，沿-Y_L方向所需控制力较大，因为在Y_L方向追踪星相对目标星有较大的初始速度，该速度由两星轨道差异造成，如果初始时刻两星位于同一轨道，则相对速度基本为0。而参考轨迹速度方向与追踪星初始速度相反，因此需要较大控制力来完成轨迹跟踪。将该控制力通过式(21)分配到追踪星本体系中，得到沿追踪星本体X_cb轴的实际控制力(图 11)。

图 12为追踪星本体X_cb轴与期望控制力的夹角。如之前所述，初始时刻所需控制力与初始速度方向相反，而追踪星初始X_cb与初始速度方向基本重合，因此，初始期望控制力与X_cb的夹角很大。通过姿态控制器，令追踪星姿态旋转，经过约50 s完成姿态跟踪。之后，追踪星X_cb指向参考轨迹的半径方向，与Z_L垂直，以此提供相应的转弯加速度。当t=8 000 s时，此时需要在Z_L方向提供推力，使追踪星沿-Z_L方向的速度减到0。因此，在Z_L方向瞬间产生了期望控制力(如图 10所示)，需要再一次进行约50°姿态机动。姿态机动完成后，X_cb指向在Z_L方向保持不变，在OX_LY_L平面匀速转动。当进入最终逼近段时，期望指向为Y_L方向，需要进行约90°姿态机动。对应的姿态角速度及控制力矩如图 13和图 14所示。

4.2 考虑初始速度方向

从图 5可以看出，如果对参考曲线的初始速度方向不加限制，则会使得追踪星初始指向与参考曲线方向差别较大，不利于轨迹跟踪。本节设计参考轨迹在初始时刻的速度方向与X_cb重合。令t₁=8 000 s，t_f可根据式(16)算得为8 086 s。期望转移轨迹及追踪星实际轨迹如图 15所示。由于将传统以Z_L为中心的螺旋线进行了旋转，使得其满足初末位置及方向要求。两次旋转角为φ=141.258°，θ=74.609°。式(14)转移轨迹参数为

$ \left\{ {\begin{array}{*{20}{l}} {A = - 6.234,B = - 6.332,r = 230.604}\\ {{D_1} = - 230.605,{D_2} = 15.599}\\ {{C_0} = - 1{\kern 1pt} {\kern 1pt} {\kern 1pt} 333.172,{C_1} = 1{\kern 1pt} {\kern 1pt} {\kern 1pt} 487.797}\\ {{C_2} = - 167.141} \end{array}} \right. $

初始时刻，参考轨线与追踪星X_cb相切，虽然不能保证速度一致，但是有利于追踪星提供切向加速度。与图 5相比，实际轨迹与参考轨迹差异更小。对比图 6与图 16，同样可以看出，当参考轨迹方向与初始X_cb重合时，跟踪误差最大为60 m，远小于任意方向时的120 m。

仿真中初始X_cb方向与初始速度方向差别很小，因此，追踪星初始速度与参考轨迹也基本相切。从图 17可以看出，速度跟踪主要是消除两者在Y_L方向的差异。对比图 18与图 10可以看出，当参考轨迹与初始X_cb相切，因为速度误差减小，轨迹跟踪所需控制力较小。除此之外，采用旋转螺旋线另一优点是整个转移轨迹只分两段，即交会段和逼近段，无需减速段，因为Z_L方向通过旋转后的螺旋轨迹可以提供相应的加速度。在8 000 s之后，追踪星到达(0, -20, 0) m的位置，进入逼近段，推力器关闭，以0.2 m/s的速度匀速驶向目标星，如图 19所示。追踪星控制力拒如图 20所示。同时，调整姿态，令X_cb指向目标，以满足对接要求，如图 21所示, 滑膜变量如图 22所示。

根据式(32)设计的姿态控制器，当姿态误差减小到一定区域时，由CLF控制切换为滑模控制。控制切换情况如图 23所示，其中1表示滑模控制，0表示CLF。CLF只作用在前100 s及后100 s，中间大部分时间为滑模控制。这种切换机制可以使得当姿态偏差较大时，以CLF保证系统最优性，节省姿态机动能量；当姿态收敛到一定范围后，切换为滑模控制，提高控制精度，增加系统对干扰的鲁棒性。

5 结论

本文研究了单方向推力器作用下航天器交会对接轨迹规划及控制问题，研究成果及结论总结如下：

1) 基于螺旋线设计了位置转移轨迹，通过引入姿态旋转矩阵，将传统螺旋线进行三维旋转，增加了转移轨迹的自由度，使其同时满足初末位置要求以及初末速度方向要求，有效降低了位置跟踪误差和位置跟踪所需控制力。

2) 为了使得追踪星本体推力方向与期望控制力方向重合，本文设计了基于CLF的姿态跟踪滑模控制策略。当姿态误差较大时，采用CLF，有效降低姿态控制能量消耗；当姿态误差较小时，切换为滑模控制，提升控制精度。

但是本文中转移轨迹设计与姿态跟踪是分开设计的，在设计转移轨迹时没有考虑期望控制力的变化速度，不利于姿态跟踪，后续有待改进。