CN1876501A - 基于行为模式的深空三轴稳定姿态定向控制方法 - Google Patents

Abstract

基于行为模式的深空三轴稳定姿态定向控制方法，涉及一种空间飞行器的姿态控制方法。为了提高星上自主能力，本发明所述方法为：带有自主导航功能的空间飞行器通过飞行计算机接收基准定向任务和辅助定向任务，空间飞行器上的飞行计算机根据导航系统提供的飞行器轨道和星历文件插值得到的目标天体轨道，确定实现上述任务组合所应满足的姿态运动学过程，将其作为参考姿态运动轨迹，采用姿态跟踪控制律，通过姿态接口单元驱动飞行器的姿态执行机构，实现需要的姿态指向。本发明的方法有利于减轻飞行器的定向负担、提高自主性。

Description

基于行为模式的深空三轴稳定姿态定向控制方法

技术领域

本发明涉及一种空间飞行器的姿态控制方法，特别是人造卫星和星际航天器的三轴姿态定向方法。

背景技术

现代深空探测器普遍采用长期三轴稳定的姿态控制方案。在科学探测期间，三轴稳定能达到很高的定向精度；此外，多数任务要求两个以上矢量的定向，三轴稳定定向是必然的解决办法。

与对地定向的地球卫星不同，单一的参考坐标系如当地轨道坐标系不能适应深空任务中多种三轴定向模式的姿态控制；此外，限于深空中尤其是小天体附近环境的多样性、未知性和复杂性，姿态敏感元件和执行机构的选配与地球卫星相比，可供选择的组合较少，当前主要依赖星敏感器、陀螺等惯性敏感器，暂时无法利用先验的天体磁场、红外等信息的测量来估计姿态。因而难以结合当前模式直接在控制回路内部产生控制指令。

因此需要一种从控制回路外部引入的参考姿态和姿态角速率的计算方法，也称姿态制导。先进、通用的星上制导算法有利于减轻飞行器的定向负担，提高自主性。Cassini、美国NEAR、欧空局Rosetta和CNES的Mars PREMIER等深空探测器都用到了姿态制导(attitude guidance)的概念。NEAR探测器依据当前的模式和飞行器、目标的位置、速度等信息，计算出即时的参考姿态四元数和角速率。在已知飞行器轨道运动的情况下，这种算法对平滑地指向运动的小天体或自旋小天体的固定特征点是非常理想的。但公开文献仅提及了最简情形下的特殊计算步骤，并没有介绍姿态运动参数的通用计算方法，从而尚没有公开的具有完整实用意义的三轴姿态控制顶层指令算法及架构。

发明内容

为了提高星上自主能力，本发明提供一种有利于减轻飞行器的定向负担、提高自主性的基于行为模式的深空三轴稳定姿态定向控制方法。

本发明按照如下方法进行深空三轴稳定姿态定向控制：带有自主导航功能的空间飞行器根据外部事件的激励，接收基准定向任务和辅助定向任务，空间飞行器上的飞行计算机根据导航系统提供的飞行器轨道和星历文件插值得到的目标天体轨道，利用下述公式确定解析参考姿态矩阵、解析参考角速度、解析参考角加速度、参考角速度、参考四元数、角速度偏差和四元数偏差，给出实现上述任务组合所应满足的姿态运动学过程，将其作为参考姿态运动轨迹，采用姿态跟踪控制律，驱动飞行器的姿态执行机构，实现需要的姿态指向，所述公式为：

解析参考姿态矩阵：C_BI＝C_BC_I ^T；

解析参考角速度： $\mathrm{\ω}_B^B{}_A=\frac{{({\stackrel{\·}{r}}_I\×{\stackrel{\·}{s}}_I)}_B}{{\stackrel{\·}{r}}_I\·s}=\frac{C_{\mathrm{BI}}({\stackrel{\·}{r}}_I\×{\stackrel{\·}{s}}_I)}{{\stackrel{\·}{r}}_I\·s_I};$

解析参考角加速度： $\stackrel{\·}{\mathrm{\ω}}_B^B{}_A=\frac{C_{\mathrm{BI}}\left\{\right[{\stackrel{\·\·}{r}}_I-2\mathrm{\ω}_B^A\×(\mathrm{\ω}_B^A\×r)]\×[{\stackrel{\·\·}{s}}_I-2\mathrm{\ω}_B^A\×(\mathrm{\ω}_B^A\×s)]{\}}_I}{{[{\stackrel{\·\·}{r}}_I-2\mathrm{\ω}_B^A\×(\mathrm{\ω}_B^A\×r)]}_I\·s_I};$

参考角速度(微分式)： ${\stackrel{\·}{\mathrm{\ω}}}_{\mathrm{ref}}=a_{\mathrm{ref}};$

参考四元数(微分式)：

角速度偏差：

四元数偏差：

姿态控制律：u＝-K_Pq_e-K_I∫q_edt-K_Dω_e。

本发明的基本原理：本发明采用导航信息和天体星历信息构造双矢量，用矢量代数的方法生成三轴姿态模式的参考姿态。参考姿态的计算类似双矢量定姿的TRIAD法，若预期两个非平行的定向矢量在本体系和参考坐标系的投影，可得参考系到本体系的坐标变换矩阵，即飞行器的参考姿态矩阵。为建立特定任务下的三轴姿态，根据任务需要，首先需给定本体系中最主要的指向矢量，称为基准轴线，对应惯性系(不单指固连于天体质心的惯性系，也可以换作动坐标系)下的定向矢量称为基准方向。基准轴线的指向由基准方向直接确定，在直角坐标系下用3个代数等式表示。指向方位与定向矢量的大小无关，相当于施加归一化约束，因此仅有2个独立等式；然后对本体的其他坐标轴作出附加的约束，选择次要定向任务的本体轴线作为辅助轴线，其绕基准轴线有一个转动自由度，对应的几何约束方程表示为1个等式。由于归一化和正交约束，姿态矩阵只有3个独立元素。由以上设定得到的方程组即可完全定义飞行器的三轴参考姿态。参考姿态角速度、角加速度通过与坐标系相对运动有关的矢量运动学运算得到。以上参考姿态参数的确定方法，为采用参考轨迹跟踪控制方法实现的空间飞行器三轴姿态控制系统的顶部指令层提供了一种基于行为模式思想的路径。该方法(或系统组件)根据外部事件的激励(一个基准定向任务和一个辅助定向任务)，给出实现该定向任务组合所应满足的姿态运动学过程，将其作为参考姿态运动轨迹，采用姿态跟踪控制律，驱动飞行器的姿态执行机构，实现需要的姿态指向。

本方法的意义在于提高星上自主能力方面，主要有以下优点：(1)主要涉及矢量线性运算，逻辑清晰，便于姿态的自主规划；(2)适合对运动状态复杂的定向目标的姿态跟踪控制；(3)根据已知的目标运动规律，直接给出参考轨迹，不需要目标的闭环跟踪；(4)在解析解的更新周期内，通过角加速度、角速度两层积分，大大降低了计算量，并解决了当两惯性矢量夹角极小情形的矢量求解奇异问题；(5)参考姿态和角速度适合采用PD、PID等简单有效的姿态控制律；必要时还可增加参考角加速度前馈；(6)提供高层控制架构，隐藏了后台操作的细节，系统行为仅取决于时间和指令，不必考虑动力学和硬件特性，便于预测；根据任务的控制精度需求、资源占用等因素，控制模式可以设计成仅在更低的执行层和硬件接口层次切换，大大减少了模式数量，简化了系统复杂程度，同时也降低了控制系统分模式设计的工作量；(7)提高了地面任务实施系统结构的各组成部分的独立性，减少相互之间的协调工作量：科学组专注于目标选择和相对于目标的指向；导航组负责提供为保持惯性矢量的数据；仪器和有效载荷组负责用最新校正数据更新本体矢量，修改当前任务的约束；(8)可用于星上的姿态运动预报，提高星务管理对几何、动力学等运动约束的监测与重构功能的自主性。

本方法相对与在先文献的创新之处在于：1、相对于本体指向轴线完全平行于本体坐标系轴线，以及要求辅助轴线与两惯性矢量共面或垂直的情形，给出了更一般意义上的姿态参考轨迹的矢量计算方法；2、给出了飞行器参考角加速度的矢量计算方法；3、利用本方法提出的解析解更新周期的概念，在该周期内通过角加速度、角速度两层积分，大大降低了参考轨迹计算的运算量，并解决了当两惯性矢量夹角极小情形的矢量求解奇异问题。

本方法的技术效果通过深空探测器多模式三轴姿态控制数学仿真与演示得到了体现。该仿真演示系统选择Ivar小行星作为探测对象，在已有的轨道设计基础上，仿真从地球逃逸轨道经转移、交会至绕飞探测的运动全过程，含一次远日点轨道机动和一次地球借力。其间的三轴姿态控制模式有对地、对日、对小行星三轴定向，以及轨道控制和自主光学导航期间的姿态机动和三轴稳定等。探测器的任务规划和调度由自主管理系统负责。当自主管理系统发送任务姿态的指令，包括目标模式和附加的定向任务，姿态控制系统就做出响应，向自主导航系统请求该时刻与目标姿态模式有关的星历服务和导航信息，计算当前的参考姿态运动参数，或者计算下一模式的三轴姿态，并执行姿态机动，转入下一模式。

在所有的三轴控制模式中，指令层都正确地给出了姿态运动参考指令，控制系统的执行层跟踪或响应该指令，使星体完成任务需要的姿态指向。

表1  典型构型的深空探测器主要任务模式下姿态制导的定向矢量

	本体基准定向轴线	空间基准定向方位	本体辅助定向轴线	空间辅助定向方位
					巡航模式	平面电池阵法线	飞行器-太阳矢量	主发动机推进轴线/中增益天线	飞行器轨道面法线/飞行器-地球矢量
自主光学导航	导航相机光轴	目标天体	平面电池阵法线	飞行器-太阳矢量
					轨道控制模式	推进系统	预先计算的推力冲量	平面电池阵法线	飞行器-太阳矢量
飞越天体模式	探测仪器视场轴线	飞行器-天体质心矢量	主发动机推进轴线/平面电池阵法线	指令方向/飞行器-太阳矢量
					绕飞小天体的科学探测模式	探测仪器视场轴线	飞行器-目标点矢量	平面电池阵法线	飞行器-太阳矢量
绕飞探测期间的对地数据传输模式	高增益天线	飞行器-地球矢量	平面电池阵法线	飞行器-太阳矢量

附图说明

图1为三轴姿态控制原理示意图，图2为姿态控制系统框图，图3为配合公式说明的单位球面矢量图，图4为本发明的整个计算与控制过程的程序流程图，图5为对日三轴定向数学仿真的偏差Euler角，图6为对日三轴定向数学仿真的角速度偏差。

具体实施方式

具体实施方式一：本实施方式需要星上导航系统和太阳系星历的支撑，根据导航系统提供的飞行器轨道确定结果和星历文件插值得到的目标天体轨道运动，包括飞行器和天体在同一时刻的位置、速度矢量。根据主、从本体定向轴线和对应的参考坐标系目标矢量，计算出从参考坐标系到本体坐标系的坐标转换矩阵，即星体的姿态矩阵。考虑相对运动造成的光行差，对定向矢量作光速修正。整个过程如图1所示。

深空探测器姿态控制系统方框图如图2所示，系统的核心是两套计算机：飞行计算机和姿态接口单元。飞行计算机作为星上姿态控制软件的载体，是姿态控制系统的中枢。从事大量的数值处理，实现多数的姿态确定、制导与控制的功能。用来实现轨道外推、敏感器信号处理、卫星姿态确定、故障判别和处理及姿态控制、轨道控制、太阳帆板指向控制和天线指向控制指令的计算等功能。姿态接口单元的主要功能是(1)提供敏感器、执行机构与飞行计算机的接口，(2)提供与星务计算机的总线联系，传递遥测、遥控及相关数据信息，(3)提供模式设置、故障记忆和健康字设置等功能。姿态接口单元搜集来自敏感器和执行器的数据传递给飞行计算机(星敏感器的测量直接输出至飞行计算机的星图识别程序)以估计当前飞行器的姿态参数，并将飞行计算机中控制器输出的控制量分配至姿态执行机构，控制姿态执行机构产生反作用力矩。此外，姿态接口单元还包含了执行安全模式的算法，故障检测、惯性测量单元重构或通过离散信号向遥测处理器报警的逻辑算法。姿态控制系统的敏感器件包括：惯性测量单元(含4个陀螺和4个加速度计，其中一套作为备份)；一个星敏感器(双探头)。姿态执行机构包括四个反作用飞轮和姿态控制推力器。惯性测量单元、星敏感器和姿态接口单元通过串口连接。姿态接口单元和飞行计算机之间通过总线交换数据。姿态控制系统利用速率陀螺和星敏感器联合确定探测器的姿态，利用反作用飞轮或者姿态控制推力器进行姿态控制，根据具体的控制模式选择使用反作用飞轮还是姿态控制推力器。姿控系统接受外部的命令和遥测处理器的姿态控制任务指令，并将当前姿控系统的状态传给命令和遥测处理器，作为遥测量，供星务管理软件分析之用。

如图4所示，本实施方式的具体姿态控制方法如下：

1、惯性矢量的管理和光行差修正

为实现姿态定向，飞行器和所有定向目标——包括天体和角尺度足够大的天体上的定点——的运动规律必须是可预报的。对带有自主功能的深空探测器来说，星载计算机主要采用圆锥曲线方程模拟太阳系天体、探测目标和飞行器的轨道运动；另外也包含了经过长期天文观测得来的太阳系大天体的精确轨道数据，用于在过渡段(圆锥曲线拼接段)采用多项式拟合提高计算准确度；通过接近飞行 阶段的自主导航，得到以探测天体为中心的坐标系下的飞行器运动规律。对于角尺度足够大的天体上的固定点的运动，需要事先了解天体的自旋矢量，这些也由自主导航系统确定。通过这些途径得到给定时刻的目标位置矢量及其矢端速度(一阶导数)。对于如行星、卫星、小天体这样的天体的质心运动加速度(位置矢量的二阶导数)，可通过简单的二体轨道力学方程近似得到：

${\stackrel{\·\·}{r}}_t=\frac{{\mathrm{\μr}}_t}{{\left|r_t\right|}^3}---\left(1\right)$

其中μ是中心天体的引力常数。

在飞行器与定向目标的轨道运动r_s/c ⁽ⁿ⁾、r_t ⁽ⁿ⁾分别已知的情况下，定向矢量及其导数由二者之差得到：

$r_{s/c-t}^{\left(n\right)}=r_t^{\left(n\right)}-r_{s/c}^{\left(n\right)},n=\mathrm{0,1,2}---\left(2\right)$

航天器，尤其是星际航天器，在参考坐标系中的相对运动速度比较大，光行差引起的角度偏差数量级最高可达分秒级。因而在需要精确定向时，须对基准定向的惯性矢量r作光速修正。记光速为c，飞行器相对于定向目标的速度是v，相对速度与飞行器—目标之间的夹角满足：

$\cos \mathrm{\θ}\frac{r\·v}{\left|r\right|\left|v\right|}---\left(3\right)$

沿飞行器相对运动方向的偏差角：

$\mathrm{\Δ\θ}={\sin }^{-1}\left(\frac{\left|v\right|\sin \mathrm{\θ}}{\left|v\right|+c\cos \mathrm{\θ}}\right)---\left(4\right)$

已知一个矢量r_a绕单位矢量n旋转σ角后得到：

r_b＝cosσ·r_a+sinσ·n×r_a+(1-cosσ)r_a·nn                            (5)

记修正前后的基准惯性矢量分别为r₀、r，利用式 (5)，则修正后的矢量：

r＝r₀cosΔθ+(n×r₀)sinΔθ                                           (6)

其中 $n=\frac{v\×r_0}{|v\×r_0|},$ 代入得：

$r=r_0\cos \mathrm{\Δ\θ}+(\frac{v\×r_0}{|v\×r_0|}\×r_0)\frac{\left|v\right|\sin \mathrm{\θ}}{\left|v\right|+c\cos \mathrm{\θ}}$

$=r_0\cos \mathrm{\Δ\θ}+(\frac{v\×r_0}{\left|v\right|\left|r_0\right|\sin \mathrm{\θ}}\×r_0)\frac{\left|v\right|\sin \mathrm{\θ}}{\left|v\right|+c\cos \mathrm{\θ}}$

$=r_0\cos \mathrm{\Δ\θ}+\frac{(v\×r_0)\×r_0}{\left|r_0\right|\left(\right|v|+c\cos \mathrm{\θ})}$

$=r_0\cos \mathrm{\Δ\θ}+\frac{(v\×r_0)\×r_0}{\left|r_0\right|\left(\right|v|+c\cos \mathrm{\θ})}$ (7)

强定向天线下行传输的角修正量与之相反： $r=r_0\cos \mathrm{\Δ\θ}-\frac{(v\×r_0)\×r_0}{\left|r_0\right|\left(\right|v|+c\cos \mathrm{\θ})}.$

假设光行差引起的矢量偏角是缓变的，为简便，在短时间内忽略其变化，对光行差造成的惯性矢量各阶导数不作修正；辅助定向矢量没有必要修正。

2、参考姿态的计算

设已知定向轴线的单位矢量r、s在飞行器本体中的坐标r_B、s_B，惯性系中对应的惯性矢量为a_I、b_I。分别在本体系和惯性系中求得两组角度关系：

$<r,s>=\arccos \left(\frac{r_B\&CenterDot;s_B}{\left|r_B\right|\left|s_B\right|}\right)=\arccos (r_B\&CenterDot;s_B)$

$<a,b>=\arccos \left(\frac{a_I\&CenterDot;b_I}{\left|a_I\right|\left|b_I\right|}\right)$ (8)

任务姿态首先要求基准轴线r对准惯性空间I中的惯性矢量a。辅助轴线s在未考虑其惯性空间的指向约束时其取值域是空间中以r＝a为旋转轴、半角为<r，s>的锥面上的单位长度的母线(图3)。显然，当s与r(即a)、b共面时，<s，b>取极值。定义下式：

φ′₊＝arccos(r·s)+arccos(a·b)

φ′_-＝|arccos(r·s)+arccos(a·b)|(9)

所有的角度限定在[0，π]：

φ₊＝min{φ′₊，2π-φ′₊}

φ_-＝min{φ′_-，2π-φ′_-}(10)

可知：

<s，b>∈[φ₁，φ₂]                                                    (11)

[φ₁，φ₂]＝[min{φ₊，φ_-}，max{φ₊，φ_-}]

假定辅助轴线s的约束可以描述为：s与b呈夹角β。β的选取取决于定向任务的需要，以适应姿态运动的内外部约束，也可为除主副定向轴线外的其他轴线的定向留出一定的弹性空间。常用的情况是β＝π/2。

列写如下的方程组求解r_I、s_I：

$r_I=\frac{a_I}{\left|a_I\right|}$

s_I·r_I＝s_B·r_B

s_I·b_I＝|s_I|b_I|cosβ＝|b_I|cosβ

|s_I|＝1                                                  (12)

该方程组的解是以r、b为旋转轴的两个圆锥的交线在坐标系I中的投影。

s_I的解有几种情况：

(1)当β∈(φ₁，φ₂)，可得两个解。根据实际情况，其中一个是伪解；

(2)仅当β取[φ₁，φ₂]的端点值时，s_I有唯一解；

(3)当β[φ₁，φ₂]，方程无解，此时可求得最接近的解。

综合起来，用下式替换式 (12)中的第三式：

s_I·b_I＝|b_I|cosφ                                        (13)

其中：

$\mathrm{\&phi;}=\left\{\begin{array}{c}\mathrm{\&beta;},\mathrm{if}{\mathrm{\&phi;}}_1<\mathrm{\&beta;}<{\mathrm{\&phi;}}_2\\ {\mathrm{\&phi;}}_1,\mathrm{if\&beta;}\&le;{\mathrm{\&phi;}}_1\\ {\mathrm{\&phi;}}_2,\mathrm{if\&beta;}\&GreaterEqual;{\mathrm{\&phi;}}_2\end{array}\right.$

s_I的求解并不通过直接求解方程组(12)(12)。首先根据几何意义，利用式(5)(5)，直接给出上面φ的后两种情况的解：

s_I(φ₁)＝cos<r·s>·r_I+sin<r·s>·(n×r_I)

s_I(φ₂)＝cos<r·s>·r_I-sin<r·s>·(n×r_I)                    (14)

其中 $n=\frac{a_I\&times;b_I}{|a_I\&times;b_I|}.$

至于一般情形s_I的解，先用球面三角形公式计算出其与s_I(φ₁)的夹角(亦等于∠srb)：

$\cos (\&angle;\mathrm{srb})=\frac{\cos \&lang;s,b\&rang;-\cos (r,s)\cos (r,b)}{\sin (r,s)\sin (r,b)}=\frac{\cos \&lang;s,b\&rang;-\cos (r,s)\cos (a,b)}{\sin (r,s)\sin (a,b)}---\left(15\right)$

然后利用式 (5)将s_I(φ₁)绕r旋转该角度得到：

s_I＝cos(∠srb)s_I(φ₁)±sin(∠srb)r_I×s_I(φ₁)+[1-cos(∠srb)]s_I(φ₁)·r_Ir_I  (16)

得到的两个解中，按下式，选取与当前测量的s_I夹角较小(即矢量内积较大)的一个作为真解。

$\cos \&lang;S_{10},S_I\&rang;=S_{10}\&CenterDot;\left[C\left(\hat{q}\right)S_I\right]---\left(17\right)$

其中 是由当前星上姿态确定给出的估计四元数

转换而来的姿态矩阵。

现在已经有了r、s分别在本体系B和惯性系I的分量，用熟知的TRIAD法可得飞行器的姿态矩阵。定义下式：

e_1B＝r_B             e_1I＝r_I

$e_{2B}=\frac{r_B\&times;s_B}{|r_B\&times;s_B|},e_{2I}=\frac{r_I\&times;s_I}{|r_I\&times;s_I|}$

$e_{3B}=\frac{e_{1B}\&times;e_{2B}}{|e_{1B}\&times;e_{2B}|},e_{3I}=\frac{e_{1I}\&times;e_{2I}}{|e_{1I}\&times;e_{2I}|}$

C_B＝[e_1Be_2Be_3B]     C_I＝[e_1Ie_2Ie_3I](18)

根据C_B＝C_BIC_I，得到本体系B相对于惯性系I的姿态矩阵：

C_BI＝C_BC_I ^T                                    (19)

矩阵C_BI可以转化为本体系B相对于惯性系I的姿态四元数。

3、参考姿态角速度的计算

对式(12)(12)求导：

${\stackrel{\&CenterDot;}{r}}_I=\frac{{\stackrel{\&CenterDot;}{a}}_I}{\left|a_I\right|}-\frac{a_I\&CenterDot;{\stackrel{\&CenterDot;}{a}}_I}{a_I\&CenterDot;a_I}{r}_I$

${\stackrel{\&CenterDot;}{s}}_I\&CenterDot;r_I+s_I\&CenterDot;{\stackrel{\&CenterDot;}{r}}_I=0$

${\stackrel{\&CenterDot;}{s}}_I\&CenterDot;b_I+s_I\&CenterDot;{\stackrel{\&CenterDot;}{b}}_I=\frac{b_I\&CenterDot;{\stackrel{\&CenterDot;}{b}}_I}{\left|b_I\right|}\cos \mathrm{\&beta;}$

${\stackrel{\&CenterDot;}{s}}_I\&CenterDot;s_I=0$ (20)

由此线性方程组解得

这里假定β是固定的。

用

表示矢量分别在两个坐标系A、B下的导数，二者之间存在如下关系

${\left({\stackrel{\&CenterDot;}{r}}_A\right)}_C={\left({\stackrel{\&CenterDot;}{r}}_B\right)}_C+{(\mathrm{\&omega;}_B^A\&times;r)}_C---\left(21\right)$

下标C表示矢量运算统一转换到某一坐标系C下进行，^Aω^B是坐标系B相对于A的角速度。

已知单位矢量r、s在本体系中是固定不变的，故：

${\stackrel{\&CenterDot;}{r}}_B={\stackrel{\&CenterDot;}{s}}_B=0---\left(22\right)$

由式(21)(21)知：

${\left({\stackrel{\&CenterDot;}{r}}_I\right)}_B={(\mathrm{\&omega;}_B^A\&times;r)}_B$

${\left({\stackrel{\&CenterDot;}{s}}_I\right)}_B={(\mathrm{\&omega;}_B^A\&times;s)}_B$ (23)

${({\stackrel{\&CenterDot;}{r}}_I\&times;{\stackrel{\&CenterDot;}{s}}_I)}_B={\left({\stackrel{\&CenterDot;}{r}}_I\right)}_B\&times;{(\mathrm{\&omega;}_B^A\&times;s)}_B$

$={\left({\stackrel{\&CenterDot;}{r}}_I\right)}_B\&CenterDot;{\left(s\right)}_B\mathrm{\&omega;}_B^B{}_A-{\left({\stackrel{\&CenterDot;}{r}}_I\right)}_B\&CenterDot;\mathrm{\&omega;}_B^B{}_A{\left(s\right)}_B$

$={({\stackrel{\&CenterDot;}{r}}_I\&CenterDot;s)}_B\mathrm{\&omega;}_B^B{}_A-{(\mathrm{\&omega;}_B^A\&times;r)}_B\&CenterDot;\mathrm{\&omega;}_B^B{}_A{\left(s\right)}_B$

$=({\stackrel{\&CenterDot;}{r}}_I\&CenterDot;s)\mathrm{\&omega;}_B^B{}_A-0$ (24)

由此解得：

$\mathrm{\&omega;}_B^B{}_A=\frac{{({\stackrel{\&CenterDot;}{r}}_I\&times;{\stackrel{\&CenterDot;}{s}}_I)}_B}{{\stackrel{\&CenterDot;}{r}}_I\&CenterDot;s}=\frac{C_{\mathrm{BI}}({\stackrel{\&CenterDot;}{r}}_I\&times;{\stackrel{\&CenterDot;}{s}}_I)}{{\stackrel{\&CenterDot;}{r}}_I\&CenterDot;s_I}---\left(25\right)$

4、参考角加速度的计算

对式(20)(20)求导：

${\stackrel{\&CenterDot;\&CenterDot;}{r}}_I=\frac{{\stackrel{\&CenterDot;\&CenterDot;}{a}}_I}{\left|a_I\right|}-2{\stackrel{\&CenterDot;}{r}}_I\frac{a_I\&CenterDot;{\stackrel{\&CenterDot;}{a}}_I}{\left|a_I\right|}-r_I\frac{{\stackrel{\&CenterDot;}{a}}_I\&CenterDot;{\stackrel{\&CenterDot;}{a}}_I+a_I\&CenterDot;{\stackrel{\&CenterDot;}{a}}_I-{(a_I\&CenterDot;{\stackrel{\&CenterDot;}{a}}_I)}^2}{{(a_I\&CenterDot;a_I)}^{\frac{3}{2}}}$

${\stackrel{\&CenterDot;\&CenterDot;}{s}}_I\&CenterDot;r_I+2{\stackrel{\&CenterDot;}{s}}_I\&CenterDot;{\stackrel{\&CenterDot;}{r}}_I+s_I\&CenterDot;{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I=0$

${\stackrel{\&CenterDot;\&CenterDot;}{s}}_I\&CenterDot;b_I+2{\stackrel{\&CenterDot;}{s}}_I\&CenterDot;{\stackrel{\&CenterDot;}{b}}_I+s_I\&CenterDot;{\stackrel{\&CenterDot;\&CenterDot;}{b}}_I=\cos \mathrm{\&beta;}[\frac{({\stackrel{\&CenterDot;}{b}}_I\&CenterDot;{\stackrel{\&CenterDot;}{b}}_I+b_I\&CenterDot;{\stackrel{\&CenterDot;\&CenterDot;}{b}}_I)}{\left|b_I\right|}-\frac{{(b_I\&CenterDot;{\stackrel{\&CenterDot;}{b}}_I)}^2}{{(b_I\&CenterDot;b_I)}^{\frac{3}{2}}}]$

${\stackrel{\&CenterDot;\&CenterDot;}{s}}_I\&CenterDot;s_I+{\stackrel{\&CenterDot;}{s}}_I\&CenterDot;{\stackrel{\&CenterDot;}{s}}_I=0$ (26)

由此线性方程组解得

由式(22)(22)、(23)(23)得：

${\left({\stackrel{\&CenterDot;\&CenterDot;}{r}}_I\right)}_B={(\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^A\&times;r+\mathrm{\&omega;}_B^A\&times;\stackrel{\&CenterDot;}{r})}_B+\mathrm{\&omega;}_B^A\&times;{(\mathrm{\&omega;}_B^A\&times;r)}_B$

$={(\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^A\&times;r)}_B+2{[\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_B$

${\left({\stackrel{\&CenterDot;\&CenterDot;}{s}}_I\right)}_B={(\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^A\&times;s)}_B+2{[\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;s)]}_B$ (27)

${[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_B\&times;{[{\stackrel{\&CenterDot;\&CenterDot;}{s}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_B$

$={[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_B\&times;{(\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^A\&times;s)}_B$

$={[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_B\&CenterDot;s_B\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^B{}_A-{[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_B\&CenterDot;\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^B{}_A{s}_B$

$={[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_I\&CenterDot;s_I\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^B{}_A-{\left(\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^A\&times;r\right)]}_B\&CenterDot;\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^B{}_A{s}_B$

$={[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_I\&CenterDot;s_I\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^B{}_A-0$ (28)

由此求得参考角加速度：

$\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^B{}_A=\frac{{[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_B\&times;{[{\stackrel{\&CenterDot;\&CenterDot;}{s}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;s)]}_B}{{[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_I\&CenterDot;s_I}$

$=\frac{C_{\mathrm{BI}}{\left\{\right[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]\&times;[{\stackrel{\&CenterDot;\&CenterDot;}{s}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;s)\left]\right\}}_I}{{[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_I\&CenterDot;s_I}$ (29)

这里用到的矢量及其一、二阶导数是不同飞行阶段下、以不同的中心天体为原点的坐标系中的矢量，由导航信息中的星历插值或目标矢量传播得到。比较方便的坐标系是原点平移到天体质心的惯性系，此时用本方法求得的仍是惯性姿态和惯性姿态角速度，不同惯性系之间的姿态之间仅是坐标变换的关系，而本体角速度不变。在惯性系和非惯性系之间，姿态需用坐标变换建立联系，而姿态角速度之间的转换还要考虑运动合成。

4、控制周期内的姿态参考值递推

前述步骤1-3用于计算一段时间内的参考姿态角速度和姿态矩阵ω_ref、C_ref(C_ref又转换为姿态四元数)的初值，这段时间称为参考值的解析计算周期，该周期是控制周期的整数倍。为减少计算量，解析计算周期内的各控制周期的参考姿态角速度通过式 (29)获得的参考角加速度a_ref的积分获得，这里假设该段时间内的角加速度是定常的；各控制周期的参考姿态四元数通过对姿态运动学微分方程的积分获得。

${\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ref}}=a_{\mathrm{ref}}---\left(30\right)$

定义 ${\mathrm{\&Omega;}}_{\mathrm{ref}}=\left[\begin{array}{c}0\\ {\mathrm{\&omega;}}_{\mathrm{ref}}\end{array}\right],$ q_ref＝[q_ref0  q_ref1  q_ref2  q_ref3]^T，

另一方面，在前述步骤1-3的解析计算中，存在一种极小的可能性，就是在某些特殊时刻，解析运算出现奇异。例如当两惯性矢量a、b之间夹角极小时，式(15)的分母趋于0。解决措施按以下的优先顺序：

(1)当参考姿态矩阵的解析计算出现奇异时，对式(31)积分，用上一个控制周期的参考角速度和姿态四元数递推算得本次参考姿态的解析参考值；

(2)当参考角速度的解析计算出现奇异时，对式(30)积分，用上一个控制周期的参考角速度递推算得本次参考角速度；

(3)当参考角加速度的解析计算出现奇异时，用上一个解析计算周期的结果代替。

而控制周期内的参考姿态参数用式(30)、(31)的数值积分方法来生成，只要给定初始条件，不会出现计算奇异的问题。因此可以避免该时刻附近更大量的计算奇异现象的出现。

设星际巡航的飞行器在本体Z轴正向装有固定的太阳帆板，X轴是主发动机推力方向。其对日定向模式的姿态定义为Z轴指向太阳，在此基础上X轴朝向推进方向。飞行器运动的原始数据如下表2、3，矢量的未知高阶导数可以通过差分方法计算。

表2

儒略日	飞行器位置矢量
				2455263.19408802	2.11083904965148	.98898784657365	.128416036168969
2455264.19408802	-2.11180350389328	.979716134928169	.129372283105007
				2455265.19408802	-2.11271754016684	.970421026515881	.130325441975227
2455266.19408802	-2.11439686139434	.951762716892242	.13222256238385
				2455267.19408802	-2.11516171672386	.942399777388388	.133166455046537
2455268.19408802	-2.11587641355616	.933014491054247	.134107189732232
				2455269.19408802	-2.11654073404841	.923606991676506	.135044731416894

表3

儒略日	飞行器速度矢量
				2455263.19408802	-5.74352083466337E-02	-.538347833669633	.055672453520749
2455264.19408802	-5.45756325896186E-02	-.539681303185491	.055497903333694
				2455265.19408802	-.051703446986874	-.541007132798628	5.53213354102787E-02
2455266.19408802	-4.59217875464246E-02	-.543636108495401	5.49622244976174E-02
				2455267.19408802	-4.30121596032857E-02	-.544939116102468	5.47796526855733E-02
2455268.19408802	-4.00898958677095E-02	-.546234376114206	5.45950504955055E-02
				2455269.19408802	-3.71549178043154E-02	-.54752181636891	5.44084030918332E-02

用前述参考轨迹的计算步骤，计算出参考姿态轨迹在各时刻的解析解如表4、5。

                                       表4

参考四元数
				0.33240536165903	-0.42889494644966	0.58633823743552	-0.60146759827216
0.33352373335111	-0.43004201553915	0.58570280213620	-0.60064799319765
				0.33464378850284	-0.43119048968220	0.58506357599421	-0.59982406459641
0.33576554468345	-0.43234038437653	0.58442052679735	-0.59899577535947
				0.33688901948797	-0.43349171510473	0.58377362195212	-0.59816308795921
0.33801423057863	-0.43464449737605	0.58312282845469	-0.59732596441312
				0.33914119573873	-0.43579874678099	0.58246811285388	-0.59648436623725

                   表5

参考角速度(×1.0e-007rad/s)	参考角加速度(×1.0e-020rad/s2)
		-0.38055528897141	-0.09154602578548
-0.38249729820893	-0.09310536545188
		-0.38447238726234	-0.09589979622418
-0.38852203842431	-0.09780847827251
		-0.39059689857002	-0.09939031319217
-0.39270531582940	-0.10097931692161
		-0.39484744230545	-0.10153832156672

姿态控制误差信号：

在参考值解析计算周期内，用角加速度及角速度的两层积分别得出参考角速度

和参考四元数q_ref。

记姿态确定过程得出的估计四元数是 星体惯性角速度在本体系的投影是

误差四元数：

式中使用了四元数乘法，上标^*代表四元数的共轭。

误差角速度：

三轴稳定姿态控制采用姿态四元数和角速度反馈的PID控制器。记误差四元数q_e＝[q_e0，q_e123]^T，取如下控制律：

u＝-K_Pq_e-K_I∫q_edt-K_Dω_e       (34)

其中q_e＝-sign(q_e0)·q_e123。

图5、图6显示了存在初始扰动和环境力矩作用下，仿真得到的姿态和姿态角速度偏差曲线。从中可见姿态跟踪达到了很高的精度，控制方法是有效的。

Claims (1)

1、基于行为模式的深空三轴稳定姿态定向控制方法，其特征在于所述方法为：带有自主导航功能的空间飞行器通过飞行计算机接收基准定向任务和辅助定向任务，空间飞行器上的飞行计算机根据导航系统提供的飞行器轨道和星历文件插值得到的目标天体轨道，利用下述公式确定实现上述任务组合所应满足的姿态运动学过程，将其作为参考姿态运动轨迹，采用姿态跟踪控制律，通过姿态接口单元驱动飞行器的姿态执行机构，实现需要的姿态指向，所述公式为：

解析参考姿态矩阵：C_BI＝C_BC_I ^T；

解析参考角加速度： $\stackrel{\&CenterDot;}{\mathrm{\&omega;}}_B^B{}_A=\frac{C_{\mathrm{BI}}{\left\{\right[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)]\&times;[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2\mathrm{\&omega;}_B^A\&times;(\mathrm{\&omega;}_B^A\&times;r)\left]\right\}}_I}{{[{\stackrel{\&CenterDot;\&CenterDot;}{r}}_I-2^A{\mathrm{\&omega;}}^B\&times;(\mathrm{\&omega;}_B^A\&times;r)]}_I\&CenterDot;s_I};$

解析参考角速度： $\mathrm{\&omega;}_B^B{}_A=\frac{{({\stackrel{\&CenterDot;}{r}}_I\&times;{\stackrel{\&CenterDot;}{s}}_I)}_B}{{\stackrel{\&CenterDot;}{r}}_I\&CenterDot;s}=\frac{C_{\mathrm{BI}}\left({\stackrel{\&CenterDot;}{r}}_I{\&times;\stackrel{\&CenterDot;}{s}}_I\right)}{{\stackrel{\&CenterDot;}{r}}_I\&CenterDot;s_I};$

参考角速度微分式： ${\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ref}}=a_{\mathrm{ref}};$

参考四元数微分式：

角速度偏差：

四元数偏差：

姿态控制律：u＝-K_Pq_e-K_I∫q_edt-K_Dω_e。

Images