基于优化的羽流撞击消旋技术瞄准策略设计

2018-07-23，，

中国空间科学技术 2018年3期

，，

中国空间技术研究院通信卫星事业部，北京 100094

在轨服务技术是空间应用技术领域的研究热点。从国内外研究现状来看，在轨服务的对象正逐渐由合作目标向非合作目标扩展。非合作目标泛指既不能提供有效合作信息，也不具备专门的在轨服务接口的航天器[1]，例如失效航天器和空间碎片。失效的GEO和MEO卫星所受的大气和地磁阻尼微弱，在太阳光压的长期作用下往往进入自由翻滚的运动状态。这类目标给在轨捕获带来很大挑战，因此国内外学者提出了多种消旋技术构想[2]，希望先为目标消旋再进行捕获。消旋技术的实现方式又可以笼统分为接触式与非接触式，非接触式消旋技术包括电磁力消旋和羽流撞击消旋[3]等。

羽流撞击消旋技术由于系统复杂度低、适用范围广、任务安全性高等优点而受到关注。羽流撞击消旋的主要构想是在服务平台部署化学推力器，令喷管指向目标卫星的太阳翼，利用羽流对目标星的撞击力形成与目标星角动量方向相反的力矩，实现对目标星的姿态调整和控制，最终为目标星消旋。相比姿态同步交会对接和机械臂抓捕等技术方案，消旋任务模式下服务平台可以避免进入目标星翻滚运动的包络区域，具有较高的安全性，且目标星上不需配置专用结构或机构，降低了系统复杂度，具有良好的适用性。

国外学者从动力学与控制角度出发开展了相关研究。Ferrari和Lavagna提出了一种羽流撞击模型和动力学仿真框架，研究了目标自旋运动情况下的消旋动力学过程，证明了技术的可实现性[4]；Thomas开展了羽流撞击消旋技术的工程探索[5]，包括服务平台构型以及轨道和姿态控制系统的设计，对化学推力器羽流模型进行了比较细致的研究，并设计了考虑服务平台轨道和姿态控制的闭环仿真系统[6]。然而上述研究存在一定的局限性，首先是研究对象均选择了ENVISAT卫星，ENVISAT的太阳翼面积较大且为矩形，利于减小羽流的非线性效应；其次ENVISAT的运动形式为绕次大惯性轴自旋，章动效应不明显，太阳翼法线与当地线速度近似重合，对羽流撞击消旋而言是很理想的工况；此外，Thomas定义瞄准策略(pointing strategy)为消旋推力器喷管相对目标应具备的位置和姿态，基于喷管轴线方向与目标的角速度矢量的几何关系给出了一种瞄准策略，该方法精度有限且受目标自身特点影响较大，不利于推广。

因此，本文重点研究以太阳翼为主要受力对象的羽流撞击消旋问题，提出了基于优化方法的羽流撞击消旋瞄准策略设计方法，重点解决一般情形下的瞄准策略。首先针对常见GEO卫星平台和太阳翼构型对消旋问题进行了归纳并给出了典型工况，其次根据参数化羽流撞击模型分析羽流撞击力相对目标位置姿态的敏感性，在此基础上提出了基于优化方法的瞄准策略规划模型，仿真结果表明该瞄准策略不但能显著提高消旋力矩精度，而且对不同的消旋工况具有良好的适用性。

1 数学模型

1.1 消旋动力学模型

高轨航天器中GEO通信卫星占据较大比例，同时由于轨位资源限制，对在轨服务的需求尤其旺盛。一般的GEO通信卫星平台采用中心体加附件的基本构型，附件包括对称或非对称安装的太阳翼和各类天线设备等。天线类设备体积较小且刚度较大，在本文中不单独考虑，简化为中心刚体的一部分。在三轴稳定卫星的动力学分析中，太阳翼一般作为柔性附件考虑，由于其展开机构处通常配置有阻尼器，因此在本文中为了降低动力学建模的复杂度，假设太阳翼为刚体，不发生弹性变形。

自由翻滚卫星通常处于整星断电状态，但太阳翼驱动机构具有一定的断电保持力矩，因此假设太阳翼保持与星本体相对转角不变，不会因羽流撞击产生的力矩转动。

在此基础上,建立卫星模型,并定义坐标系如下：

{b},目标本体坐标系，原点位于星体质心，zb沿星体纵向，xb轴为太阳翼安装方向，yb轴由右手坐标系确定。

{A},附件坐标系，初始状态下，原点位于星体质心，xA轴为太阳翼的纵向，zA轴沿太阳翼法线方向且与星本体zb轴重合，yA轴沿太阳翼横向。

单轴式太阳翼可绕xA轴转动，在初始安装位置下，xA轴与xb轴重合。目标星基本构型以及坐标系关系如图1所示。

图1 目标星构型Fig.1 Configuration of target

根据最大惯性轴定理，非对称刚体的自由运动趋于演变为绕最大惯量轴的自旋运动[7]。一些针对在轨失效航天器的研究也表明，失效卫星往往形成绕最大或者次大主惯性轴的近似自旋运动[6]。基于以上分析，针对本文选取的卫星构型，假设卫星近似自旋运动，角动量矢量指向星本体xb轴附近。刚体卫星的姿态动力学方程为:

(1)

式中：ω为卫星瞬时角速度；I为目标星惯量阵；T为绕质心的力矩。假定目标卫星本体系与惯性主轴重合，则在{b}系下I为对角矩阵。

根据文献[6]，期望消旋力矩矢量可以按照如下方法确定：依据目标星姿态运动参数，选取适当的控制器来确定期望角加速度，从而根据式(1)得到期望消旋力矩矢量。针对于处于自旋状态的目标星，考虑比例控制则期望角加速度为：

(2)

式中：Kp为比例系数；ωref为控制目标，在消旋问题中应令ωref=0。考虑太阳翼相对初始位置的转角a，{b}到{A}的转换矩阵为：

(3)

显然，在{A}中研究消旋推力器相对太阳翼的位置姿态关系更加直观简便，可以直接建立期望消旋力矩与消旋推力器喷管理想位置姿态之间的对应关系，因此下文中如不加特殊说明，标量运算均在{A}中进行。

为了简化问题便于分析计算，本文忽略太阳翼的表面特征，认为太阳翼为连续的无厚度矩形平板。在消旋操作过程中，总是对某一侧太阳翼施加羽流撞击作用，因此不妨以+xA一侧的太阳翼为研究对象，尺寸参数如表1所示。

表1 太阳翼几何参数

1.2 羽流场模型

Simons模型[8]是一种参数化羽流分布模型。模型指出，羽流传播距离远大于喷管出口截面尺寸时，羽流场中任一点处的压强可以用流场的几何参数表示。基于分析建立的羽流密度分布函数为：

(4)

式中：ρ*为喷管喉部的压强；ρ为流场中任一点处的羽流密度大小；r为该点距离喷管出口的距离标量；R*为喷管喉部半径；f(θ)为羽流角分布函数，θ为沿喷管出口到该点视线方向与推力器轴线方向的夹角；羽流常数Ap由气体性质所决定，

(5)

式中：θlim为推力器羽流的极限偏转角；U*/Ulim为气体极限膨胀状态下的当地马赫数。根据等熵流动特性，得到近似关系：

(6)

式中：γ为气体比热比。

羽流的角分布函数f(θ)为：

(7)

式中：θ∞为超音速无粘流时的θlim值。Simons考虑喷管出口处附面层效应得到羽流随偏转角度增大的衰减满足指数规律[7]：

f(θ)=f(θ0)e-β(θ-θ0)

(8)

式中：θ0和β的取值取决于喷管出口处的状态。如果流动是理想的无粘流，附面层中的流动也符合余弦规律：

(9)

式中：δ为喷管出口处的附面层厚度；Re为喷管出口处的截面半径。当θ0接近θ∞时，有以下近似关系成立：

(10)

由此可得β：

(11)

1.3 表面作用模型

羽流充分膨胀后形成自由分子流动，在与物面接触的过程中通过碰撞、反弹和吸附等行为将动量传递给目标。Feshe模型对于羽流动量传播采用的假设为：沿喷管出口轴线方向上微元受力应该等于喷管出口处的压力。文献[6]对羽流模型进行了分析比较，将Simons模型的分布规律代入Feshe模型：

(12)

式中：C为常数，

(13)

羽流气体与物体表面的相互作用是复杂的物理过程，通常用镜面反射系数和散射系数来描述传播的动量[9]：

(14)

式中：Uw为完全热适应情况下分子离开壁面时的反射速度；v为入射角；n为当地法向；ev为来流方向;σn和στ分别为法向和切向的动量传播系数，根据物体表面特性和羽流特性确定。

本文中，根据LEO轨道上试验数据[10]给出以下假设：

1)化学推力器的羽流是高热自由分子流动；

2)壁面温度300 K，得到速度比为0.23；

3)气体采用完全热适应假设，存在97%漫反射。

针对消旋推力器参数的选取，本文以某单组元10 N推力器作为参考，给出设计参数如表2所示。

表2 消旋推力器设计参数

2 消旋策略

2.1 敏感性分析方法

除了在轨服务技术中普遍存在的技术难点[11-13]外，羽流撞击消旋技术的复杂性主要体现在羽流模型的非线性。在消旋问题中这种非线性体现在两个方面：第一，随着入射流的偏转以及目标表面的特征，羽流撞击所产生合力的大小和方向存在非线性；第二，由于羽流的扩散效应以及目标尺寸有限，导致总有部分羽流不能被目标捕获而造成非线性。羽流敏感性分析是研究消旋力矩相对喷管相对位置姿态的变化趋势，对非线性进行初步分析，从而为瞄准策略设计提供基础。

根据式(14)给出的流场中某一面积微元所受羽流撞击力，在附件系下对太阳翼表面区域积分得到太阳翼所受的羽流撞击力为：

(15)

表面力作用的积分属于第一类曲面积分，可以转化为每个分量在该区域上的二重积分，即：

(16)

类似有羽流撞击力所产生的相对目标星质心的力矩为:

(17)

因此消旋力矩可以表示为喷管相对位置姿态的函数：

T=T(xN,yN,zN,α,β)

(18)

式中：(xN,yN,zN)为消旋推力器喷管的坐标；(α,β)为推力器喷管指向的方位角，定义方式与极坐标相同，偏航角α∈[0,2π]表示喷管指向矢量与xA轴正向夹角，俯仰角β∈[0,π]表示喷管指向矢量与zA轴正向夹角。

由式(18)可知消旋力矩矢量由5个变量确定，对应喷管相对太阳翼的5自由度运动。因此在敏感性分析中分别考虑推力器典型的平动和转动对消旋力矩的影响。

平动情况下令喷管垂直指向太阳翼，即β=π，喷管沿3个坐标轴平动，从而分析传播距离和捕获面积对消旋力矩的影响。转动情况下，选取太阳翼中心区域某一定点作为瞄准点，改变喷管轴线相对指向，从而研究在不同入射方位下的消旋力矩变化规律。

2.2 瞄准优化方法

(1)约束条件

在整个消旋任务中，服务星应该满足安全性约束，本文的安全约束区域按照球形安全区域设计，给出极限安全距离约束条件为:

(19)

式中：RN为在附件坐标系下的位置矢量；dmin为最小安全半径。

(2)目标函数

根据分析，瞄准策略的设计可用优化问题描述，优化目标是使消旋力矩的方向尽量接近期望方向的同时使幅值最大化。根据指数函数特点不妨令目标函数为：

(20)

式中：B为常数;‖T‖表示消旋力矩矢量的范数，

(21)

ε为实际力矩与期望力矩方向的夹角，即误差角:

(22)

该目标函数具备以下性质：

1)消旋力矩幅值‖T‖增大以及误差夹角ε均使得目标函数降低。

2)在‖T‖增大的初段，目标函数下降迅速，此时力矩幅值所占的权重较大；在‖T‖增长的末段，指数项和三角函数项的变化速度相近，此时精度与幅值具有相近的权重。

3)调控参数B可以改善函数特性，对优化结果有一定影响。

结合约束条件得到瞄准策略的优化模型为:

(23)

3 仿真分析

3.1 敏感性仿真

根据第2.1节中的分析，可以设计仿真如下：仅平动情况为喷管轴线始终指向-z轴方向，控制变量β=π。消旋力矩相对喷管坐标的变化规律如图2所示。

图2(a)为喷管坐标y=0时，不同z向距离下，Ty相对x坐标的变化规律，随着羽流传播距离越短，Ty极值点对应的x坐标逐渐增大。由于太阳翼关于x轴的对称性，Tx和Tz的大小不受x改变影响。同时，Ty最大值与z的大小近似成平方反比。

图2(b)(d)分别展示了在不同z向距离下，Tx和Tz相对y坐标的变化规律。其中Tx和Tz均随y增大而增大，羽流撞击的切向力所产生的力矩要大于等效力臂增大产生的力矩，从而说明针对当前构型的太阳翼，羽流撞击消旋能够很好地适应太阳翼倾斜的工况。

图2(c)表示Ty相对y坐标的变化规律。Ty随着y增大而减小，可以认为是由于太阳翼捕获的羽流减少了。

此外，可以看出Tx量级小于Ty和Tz，这是由于当喷管指向位置沿+y偏移时，切向合力产生+x方向的力矩，法向合力产生-x方向的力矩，两者存在一定的抵消效应，这种效应使得绕太阳翼x轴方向的转动相比绕另外两轴的转动运动更难消除。另一方面也可以认为羽流撞击的方式产生的Tx较小，能够满足太阳翼不发生转动的基本假设。

其次分析消旋力矩相对入射角的变化规律。令喷管始终指向点(9.8,0,1.5)，传播距离为8 m，在不同的喷管指向角下，根据喷管相对位置姿态关系，可以反解出喷管的坐标，从而可以分别计算喷管沿不同指向时消旋力矩各分量的变化规律。利用对称性可以降低计算量，当β∈[90°,180°]时，表示喷管垂直向下，力矩变化规律与β∈[0,90°]时对称；同理，在a∈[0°,180°]时的力矩分布规律与a∈[180°,270°]对称，因此只考虑a∈[0°,180°]，β∈[120°,180°]时力矩变化规律。

图3展示了3个方向的力矩分量相对喷管方位角的变化规律。从仿真结果可以看出，推力器越接近于垂直方向，即俯仰角β接近180°时，Ty分量到达最大，同时合力矩最大。当喷管在yoz平面附近以较大入射角施加羽流时，偏航角α在90°或270°附近，而俯仰角β较小，此时对应的羽流捕获面积最小的情况，合力矩明显减小，由于y方向的羽流撞击力增大，此时Tz有最大值，意味着改变喷管倾角能够产生相当的切向力从而产生较大的z向消旋力矩。此外，在偏转入射时Tx的量级同样小于Ty和Tz，说明倾斜入射工况下Tx的抵消作用仍旧存在，与之前的分析一致。

图2 平动时消旋力矩变化规律Fig.2 Torque relates to translational motions

图3 喷管转动对消旋力矩影响Fig.3 Influence of thruster rotation on torques

综合以上分析可以得到如下结论：羽流撞击所产生的消旋力矩能够产生足够大的Ty分量，是消旋实现的核心要求；喷管倾斜时能够产生同一量级的Ty和Tz，较难产生相同量级的Tx，这一特点符合羽流撞击消旋的任务需求，同时对于任务规划具有指导作用。在羽流撞击消旋的过程中，应当尽量保持喷管轴线指向目标平板的横向中心区域以减轻非线性效应。

3.2 消旋仿真

由于在实际消旋过程中，化学推力器羽流的传播速度远大于服务平台与目标之间的相对运动速度，因此忽略羽流场的瞬态效应以及双星相对运动造成的相对位置姿态改变。同时，羽流撞击消旋总是针对一侧的太阳翼开展的。随着目标周期性转动运动，总以逐渐接近的太阳帆板为主要受力对象，因此假设目标其他部分受到的羽流撞击忽略不计。

根据式(23)所示的瞄准策略规划模型，设计仿真算例如表3所示。情况1为期望消旋力矩方向Tdes在y轴附近，情况2和3为z方向期望消旋力矩分量较大；同时为考察目标函数参数不同所造成的影响，令B取值为0.01和0.5。

表3 仿真参数设置

仿真计算的结果如表4和图4所示。对于期望力矩方向接近y轴正向的情况，该优化方法能够使实际消旋力矩与期望方向夹角收敛到0°，同时能保持较大的力矩幅值，这对于消旋是非常有利的。对于Tz较大的情况，调整喷管相对位置姿态能够产生足够的z方向消旋力矩。此时可能出现问题是，力矩的幅值和指向精度需要权衡，如图4(b)所示，追求较小的误差角度可能导致力矩幅值减小，此时可以通过合理配置参数A进行修正。

从仿真结果中可以看出A对精度和消旋力矩幅值有明显的调整作用，因此有必要分析参数A的影响。保持期望消旋力矩矢量为(-0.01,1,0.4)，令A在[0.01,10]的范围内变化，得到参数A对消旋力矩各分量以及模值大小Tabc的变化规律如图5所示。参数A的取值对消旋力矩有明显影响，随着A取值的减小，消旋力矩的方向误差和消旋力矩的幅值均减小，相当于增加了角度误差的权重，使得寻优结果倾向于角度误差更小。当A逐步增大，当力矩幅值收敛时，角度误差约为12°。因此，可以认为A的理想取值范围应大约在[1,10]，在保证力矩较大的情况下有足够精度。