基于大气阻力的微小卫星编队飞行控制方法及控制装置

摘要

基于大气阻力的微小卫星编队飞行控制方法及控制装置，涉及微型航天器与航天测控技术领域。本发明的GPS接收机的信号输出端与轨道确定装置的信号输入端连接，太阳电池阵列和三轴磁强计的信号输出端与姿态融合滤波的输入端连接，三轴磁强计的输出端与纯磁测定姿的输入端连接，姿态融合滤波的信号输出端与星间链路的第二输入端连接，星间链路的第一输出端与编队队形控制律的输入端连接，星间链路的第二输出端与卫星姿态控制律的输入端连接，通过控制气动板来进行编队卫星飞行的队形和姿态控制。本发明实现了延长了航天器的寿命，功耗低、成本低、无污染的基于大气阻力的目的。

说明书

技术领域

本发明涉及微型航天器与航天测控技术领域。

背景技术

微小卫星技术近十年来一直是国内外航天领域研究的热点，而与此相关的若干关键技术的研究已经吸引了包括高校在内的众多研究机构的关注，由于微小卫星具有重量轻、体积小、成本低、研制周期短以及功能密度高等一系列优点，在通讯、遥感、军事、行星探测、工程技术实验等领域发挥着重要作用。研究开发具有高指向精度和稳定度的姿态测控技术具有十分重要的意义。

编队飞行是微小卫星的一个重要发展方向，其兴起标志着航天技术进入一个崭新的时代，微小卫星编队飞行是由多颗彼此之间距离较近的卫星组成，具有特定的空间几何构形，通过星间链路实现协同工作，不仅能完成特定的航天任务，而且在体积、性能、费用等方面比单颗大卫星具有明显的优势。

为了较好的实现卫星编队功能，编队飞行的各卫星必须密切配合协同工作，其姿态和轨道控制是实现编队飞行的关键技术之一。传统的卫星姿态控制通常采用动量轮或推进系统提供控制力矩，编队飞行卫星的队形控制通常使用推进系统来提供推力。现有的技术主要存在以下三个问题：一是采用动量轮或者推进系统的控制系统不仅技术复杂而且成本；二是微小卫星携带燃料有限，姿态和队形控制需要燃料的持续消耗，一旦耗尽卫星就失去了价值；三是推进器的羽流对卫星的冲击不仅使其稳定性受到了干扰，而且对精密的星载仪器造成污染。本发明给出了一种利用大气阻力进行微小卫星编队飞行控制的装置，给出了相应的系统配置和姿态轨道控制方法，使用本发明将使得卫星的编队飞行任务将不受限于卫星所载燃料的多少，延长了航天器的寿命，尤其适用于低功耗、低成本、轻重量、无污染的卫星编队飞行任务，具有很好的应用前景。

发明内容

本发明目的是提供一种利用大气阻力进行微小卫星编队飞行轨道和姿态控制的方法，延长了航天器的寿命，功耗低、成本低、无污染的基于大气阻力的微小卫星编队飞行控制方法及控制装置。

本发明为实现上述目的，采用如下技术方案：

本发明包括GPS接收机、轨道确定装置、星间链路、太阳电池阵列、姿态融合滤波器、三轴磁强计、纯磁测定姿、编队队形控制律、卫星姿态控制律、气动板；GPS接收机的信号输出端与轨道确定装置的信号输入端连接，将GPS接收机接收的轨道信号经过轨道确定装置计算后得到准确的轨道状态；太阳电池阵列和三轴磁强计的信号输出端与姿态融合滤波的输入端连接，姿态融合滤波器对太阳电池阵列测量的太阳方位信息和磁强计测量的地球磁场信息进行滤波，得到高精度姿态信息；三轴磁强计的输出端与纯磁测定姿的输入端连接，将三轴磁强计采集的磁场信息经过纯磁测定姿的计算后得到卫星姿态信息；姿态融合滤波的信号输出端与星间链路的第二输入端连接，纯磁测定姿态的信号输出端与星间链路的第三输入端连接，经过计算得到各卫星相对和绝对的轨道与姿态信息；星间链路的第一输出端与编队队形控制律的输入端连接，将星间链路中的准确的轨道信息传输到编队队形律中；星间链路的第二输出端与卫星姿态控制律的输入端连接，将星间链路中的准确的姿态信息传输到卫星姿态控制律中；编队队形控制律、卫星姿态控制律的输出端分别与气动板连接，通过控制气动板来进行编队卫星飞行的队形和姿态控制。

比较好的是，本发明的小卫星为立方体结构，气动板设置在小卫星飞行方向的背面，气动板与小卫星的飞行方向垂直，通过控制气动板的伸缩来改变卫星的迎风面，一方面产生气动力矩来控制卫星姿态，另一方面使得编队飞行的各卫星之间产生差动大气阻力，用来对编队卫星的队形进行控制。

本发明基于大气阻力的微小卫星编队飞行的控制装置的控制方法：

第一步：卫星轨道确定。通过GPS接收机的信号确定卫星的轨道信息；若GPS信号丢失，则采用广义卡尔曼滤波方法进行轨道信息的递推；

第二步：卫星姿态确定；

第三步：编队飞行卫星相对轨道和姿态确定：经过第一步和第二步得到各颗卫星的轨道和姿态信息后，通过星间链路得到各卫星间的相对和绝对的轨道和姿态信息；

第四步：编队卫星姿态控制：上述第三步的相对和绝对姿态信息传输到卫星姿态控制律中，卫星姿态控制律控制气动板与大气的有效接触面积，应用气动力矩进行滚动、俯仰和偏航方向的姿态；

第五步：编队飞行队形控制：上述第三步的相对和绝对轨道信息传输到编队队形控制律中，并判断第四步的的姿态控制是否完成，若姿态控制完成后通过卫星编队队形控制律驱动气动板与大气的有效接触面积，应用各卫星间的大气差动阻力进行编队卫星的队形控制；

第六步：星间链路对各卫星之间的姿态信息与设定的姿态参数进行比对判断，若判断完成姿态控制任务，则继续监测；若判断为否，则执行上述第四步；

第七步：星间链路对各卫星之间的轨道信息与设定的编队飞行的队形进行判断，若判断完成队形控制任务，则继续监测；若判断为否，则执行上述第五步。

比较好的是，本发明的卫星姿态确定存在的以下三种情况：

情况1：当仅仅需要角速率信息时，利用三轴磁强计测量的磁场信息，通过纯磁测定，采用磁场差分的B-dot方法来间接的确定卫星的角速率变化；

情况2：当同时需要卫星的三轴姿态和角速率信息，并且卫星处于太阳光照区时，利用太阳电池阵列、三轴磁强计测量的太阳方位和地球磁场信息，以及经过第一步得到的卫星轨道信息，通过姿态融合滤波器来确定高精度的卫星姿态；

情况3：当同时需要卫星的三轴姿态和角速率信息，并且卫星处于地球阴影区时，利用三轴磁强计测量的磁场信息，以及经过第一步得到的卫星轨道信息，通过纯磁测定姿，采用EKF的方法来确定卫星的姿态信息。

比较好的是，本发明的第二步的卫星姿态确定方法，采用太阳电池阵列和磁强计组合，通过卫星动力学方程计算角速度来设计无陀螺滤波器，滤波器采用双矢量定姿q-method方法和非线性EKF滤波相融合，首先应用q-method方法将太阳矢量和磁场矢量的六维观测量转换为四元数，降低观测维数并将观测方程转换为线性，然后代入EKF进行滤波计算，得到卫星的姿态信息。

比较好的是，本发明第四步的编队飞行卫星的姿态控制，仅仅采用气动板作为执行机构，根据纯磁测差分计算的磁场变化率信息设计卫星速率阻尼控制律，根据姿态确定算法得到的姿态四元数信息和角速度信息设计卫星姿态捕获及稳定控制律，得到期望的气动力矩后进一步设计气动板驱动控制律，从而控制编队飞行卫星的姿态。

比较好的是，本发明第五步编队飞行卫星的队形控制仅仅采用大气阻力，根据编队飞行的各卫星的相对位置信息，通过建立编队卫星的相对运动模型，针对共面跟飞和共面椭圆编队的终端约束，设计最优控制律进行卫星编队飞行的队形建立、队形重构以及队形保持等功能控制。

本发明的基于大气阻力的微小卫星编队飞行控制系统，在卫星入轨后的主要工作模式及工作流程为：

1、消除星箭分离引起的初始扰动

对于编队飞行的各卫星，利用卫星本体所处的磁场变化情况，采用磁强计测量的磁场变化信息，间接的测量卫星角速率变化，采用设计的磁场变化率反馈控制律来控制气动板，应用气动力矩进行滚动、俯仰和偏航方向的速率阻尼。

2、姿态捕获及三轴姿态稳定控制

对于编队飞行的各卫星，在速率阻尼完成后，光照区时启动太阳电池阵列加磁强计的姿态确定组合滤波器，阴影区时采用纯磁测定姿方法，得到各颗卫星的绝对姿态后可以根据星间链路计算卫星间的相对姿态，根据任务的需要采用相应的姿态反馈控制律驱动气动板，应用气动力矩进行卫星各轴的绝对及相对的姿态控制。

3、卫星编队飞行的队形建立控制

对于编队飞行的各卫星建立姿态基准后，根据编队队形的任务要求，应用大气阻力建立编队队形。由于大气阻力仅可以提供沿航迹向的作用力，而不能提供轨道面法向的作用力，因此本发明的基于大气阻力编队队形控制主要针对共面编队飞行的情况，如共面的跟飞、绕飞、椭圆等编队控制。根据初始条件和最终目标设计最优控制律，根据卫星间的差动大气阻力形成期望的编队队形。

4、卫星编队飞行的队形保持控制

编队飞行的卫星，在轨运行期间会受到各种摄动力，使之不能保持预设的运行轨迹，从而破坏整个编队队形。根据各卫星间的相对轨道和姿态信息，设计最优控制律，通过控制气动板来保持编队卫星间的相对位置，从而保证编队队形的稳定性。

5、卫星编队飞行的队形重构控制

对于编队飞行的各卫星，根据任务的要求需要进行队形重构，例如对于共面编队飞行，根据需要控制卫星跟飞距离的增大或减小，控制绕飞椭圆的膨胀与收缩等，根据任务要求设计最优控制律，通过控制气动板来进行共面编队队形的重构。

本发明采用上述技术方案，与现有技术相比具有如下优点：

1.国内外卫星的轨道控制通常都采用推进系统，不仅系统复杂，而且需要燃料消耗。微小卫星由于体积、质量、功耗等限制，携带燃料有限，一旦耗尽卫星就失去了价值。同时推进器的羽流不仅会对星载仪器造成污染，而且对卫星造成冲击使其稳定性受到干扰。采用本发明可以使编队卫星飞行时间不受限于所载燃料的多少，尤其适用于低功耗、低成本、轻重量、无污染的卫星编队飞行任务。

2.对于编队飞行的相对和绝对姿态控制，以及单颗卫星的姿态控制，通常的姿态执行机构采用推进系统、动量轮(包括反作用轮和偏置动量轮)以及磁力矩器等，或者采用不同的执行机构组合，姿态控制系统设计较为复杂。本发明的姿态控制仅仅使用气动板，不仅是一种新的姿态控制方法，而且姿态控制的同时还可以进行编队控制，一物多用提高功能密度，显著减小控制系统的复杂性，符合低成本、低功耗的微小卫星设计思想。

3.本发明设计了一种融合姿态确定性算法q-method和非线性最优滤波EKF的双重姿态确定滤波器，根据太阳电池阵列和磁强计测量，结合磁场模型和太阳矢量模型进行姿态确定，不仅大大减小采用直接观测方法的计算量，而且可以提供姿态确定精度。本发明设计了基于大气阻力的编队卫星姿态控制和队形控制律，有效的进行微小卫星编队飞行控制。本发明中的姿态确定方法与编队控制方法，不仅应用于基于大气阻力的编队飞行控制，而且可以为采用其它控制方式的编队卫星控制也提供了重要参考。

4.本发明的大气阻力进行微小卫星编队飞行的控制系统，具有成本低、重量轻、功耗小、无污染、寿命长的优点，完全符合微小卫星的“快、好、省”的发展原则。本发明中采用的各项关键技术及创新思想，可以推广到空间飞行器的轨道和姿态控制系统的设计中，在中、低轨道卫星编队飞行、星座系统控制以及单颗卫星姿态轨道控制等领域的研究中，具有非常广阔的工程应用前景。

附图说明

图1是本发明的基本原理图。

图2a是本发明的卫星结构示意图。

图2b是图2的俯视图。

图2c是图2的左视图。

图2d是图2的侧视图。

图3是本发明的飞行姿态确定滤波器流程。

图4a是本发明的共面跟飞编队状态示意图。

图4b是本发明的共面椭圆编队状态示意图。

图5a是本发明的卫星飞行姿态确定仿真结果的参考运动轨迹示意图。

图5b是本发明的卫星飞行姿态确定仿真结果的直接观测EKF时四元数误差示意图。

图5c是本发明的卫星飞行姿态确定仿真结果的双矢量定姿时四元数误差示意图。

图5d是本发明的卫星飞行姿态确定仿真结果的双重滤波器时四元数误差示意图。

图5e是本发明的卫星飞行姿态控制仿真结果的速率阻尼控制角速度结果示意图。

图5f是本发明的卫星飞行姿态控制仿真结果的速率阻尼控制时气动板长度变化示意图。

图5g是本发明的卫星飞行姿态控制仿真结果的三轴姿态控制四元数结果示意图。

图5h是本发明的卫星飞行姿态控制仿真结果的三轴姿态控制时气动板长度变化示意图。

图5i是本发明的卫星飞行队形控制仿真结果的共面跟飞编队队形建立示意图。

图5j是本发明的卫星飞行队形控制仿真结果的共面椭圆编队队形建立示意图。

图5k是本发明的卫星飞行队形控制仿真结果的共面跟飞→椭圆编队重构示意图。

图5l是本发明的卫星飞行队形控制仿真结果的共面椭圆→跟飞编队重构示意图。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

1、技术实施整体路线：

如图1所示，本发明包括GPS接收机11、轨道确定装置12、星间链路13、太阳电池阵列14、姿态融合滤波器15、三轴磁强计16、纯磁测定姿17、编队队形律18、卫星姿态控制律19、气动板；GPS接收机11的信号输出端与轨道确定装置12的信号输入端连接，将GPS接收机11的信号经过轨道确定装置12计算后得到准确的轨道状态，轨道确定装置12的第一信号输出端与星间链路13的第一信号输入端连接；太阳电池阵列14的信号输出端与姿态融合滤波器15的第一信号输入端连接，三轴磁强计16的信号输出端与姿态融合滤波器15的第二信号输入端连接，姿态融合滤波器15对太阳电池阵列14测量的太阳方位信息以及三轴磁强计16测量的磁场信息进行融合滤波计算后得到高精度姿态信息；三轴磁强计16的输出端与纯磁测定姿17的输入端连接，将三轴磁强计16测量的磁场信息经纯磁测定姿17的计算后得到姿态信息；姿态融合滤波器15的信号输出端与星间链路13的第二输入端连接，纯磁测定姿17的信号输出端与星间链路13的第三输入端连接；星间链路13的第一输出端与编队队形控制律18的输入端连接，将经过星间链路13处理后得到的相对和绝对的轨道信息传输到编队队形控制律18中；星间链路13的第二输出端与卫星姿态控制律19的输入端连接，将星间链路13处理后得到的相对和绝对姿态信息传输到卫星姿态控制律19中；编队队形控制律18、卫星姿态控制律19的输出端分别与气动板连接，通过控制气动板来进行编队卫星飞行的队形和姿态控制。

2、卫星设计、大气阻尼及气动力矩分析：

上述基于大气阻力的微小卫星编队飞行控制系统所设计的卫星结构示意如附图2a、2b、2c、2d所示。与目前国际上大多数的小卫星、纳卫星以及皮卫星类似，卫星采用长方体结构，所不同的是在飞行方向的背面上安装四块轻型的可以伸缩的气动板，作为姿态控制和队形控制的执行机构。记卫星长、宽和高的尺寸分别为a、b和c，气动板宽为e，通过调节长度d来控制气动板的面积。对于姿态控制坐标系定义为：轨道系F_o原点位于卫星质心，飞行方向为x轴，y轴指向轨道面负法向，z轴按右手规则定义；星体系F_b原点位于卫星质心，三轴姿态稳定时与轨道系重合；惯性系F_i原点位于地心，x轴指向春分点z轴沿地球自转轴，y轴由右手规则定义。给卫星各个外表面进行编号：1和2分别对应x方向的两个面，3和4对应y轴，5和6对应z轴，4块气动板分别编号为7、8、9和10。则卫星六个面的中心以及四块气动板的中心相对于卫星质心的矢径分别为：

${\overrightarrow{r}}_1=\left(\begin{array}{c}a/2\\ 0\\ 0\end{array}\right),$${\overrightarrow{r}}_2=\left(\begin{array}{c}-a/2\\ 0\\ 0\end{array}\right),$${\overrightarrow{r}}_3=\left(\begin{array}{c}0\\ b/2\\ 0\end{array}\right),$${\overrightarrow{r}}_4=\left(\begin{array}{c}0\\ -b/2\\ 0\end{array}\right),$${\overrightarrow{r}}_5=\left(\begin{array}{c}0\\ 0\\ c/2\end{array}\right),$${\overrightarrow{r}}_6=\left(\begin{array}{c}0\\ 0\\ -c/2\end{array}\right)---\left(1\right)$

${\overrightarrow{r}}_7=\left(\begin{array}{c}-a/2\\ e/2\\ -(c+d_1)/2\end{array}\right),$${\overrightarrow{r}}_8=\left(\begin{array}{c}-a/2\\ -e/2\\ -(c+d_2)/2\end{array}\right),$${\overrightarrow{r}}_9=\left(\begin{array}{c}-a/2\\ e/2\\ (c+d_1)/2\end{array}\right),$${\overrightarrow{r}}_{10}=\left(\begin{array}{c}-a/2\\ -e/2\\ (c+d_1)/2\end{array}\right)---\left(2\right)$

当卫星轨道高度在120km以上时，大气运动可以看作自由分子流，此时相对于轨道系的大气阻力和气动力矩公式分别为

$\overrightarrow{F}=c_d/2\left(\rho V_R^{2}A\right)\overrightarrow{v},$$T=1/2c_d\rho V_R^{2}A(\overrightarrow{r}\times \overrightarrow{v})---\left(3\right)$

其中，c_d为阻力系数，ρ为大气密度，A为迎风面面积，为来流的单位矢量，为气动中心到卫星质心的矢径，对于本发明所设计的卫星外形10个面的矢径分别见式(1)、(2)。V_R为卫星与大气之间的相对速度，即

$V_R^2=\frac{\mu }{R_c}(1-\frac{3{\omega }_e}{{\omega }_o}\cos i)---\left(4\right)$

其中，μ为地球引力常数，R_c为卫星到地球中心的距离，ω_e为地球自转角速度，ω_o和i分别为轨道角速度和轨道倾角。来流方向单位矢量在轨道坐标系分量为

$\overrightarrow{v}={\left(\begin{array}{ccc}-1& a& 0\end{array}\right)}^T---\left(5\right)$

其中，a＝(1.5ω_e/ω_o)sin i cos ω_ot，t为从升交点开始计算的时间。大气密度计算采用指数衰减大气模型

ρ＝ρ_o exp(-(h-h_o)/H)                                 (6)

其中，ρ_o为参考高度为h_o时的参考大气密度值，h为实际高度，H为标高。

对于本发明的编队卫星，对卫星每一个面所受的大气阻尼和气动力矩进行计算，然后进行求和得到整个卫星所受大气阻尼和气动力矩分别为

${\overrightarrow{F}}_a=\underset{n=1}{\overset{6}{\Sigma }}{\overrightarrow{F}}_n+\underset{m=7}{\overset{10}{\Sigma }}{\overrightarrow{F}}_m,$$T_a^o=\underset{n=1}{\overset{6}{\Sigma }}{\overrightarrow{T}}_n+\underset{m=7}{\overset{10}{\Sigma }}{\overrightarrow{T}}_m---\left(7\right)$

3、编队卫星的轨道姿态确定方法设计：

本发明的微小卫星编队飞行姿态确定滤波器设计流程如附图3所示。采用四元数$q={\left(\begin{array}{cc}{q}_0& {\overrightarrow{q}}^T\end{array}\right)}^T$描述卫星的姿态，则轨道系到星体系的姿态旋转矩阵为

$C_o^b=\left(\begin{array}{ccc}1-2(q_2^2+q_3^2)& 2(q_1{q}_2+q_0{q}_3)& 2(q_1{q}_3-q_0{q}_2)\\ 2(q_1{q}_2-q_0{q}_3)& 1-2(q_1^2+q_3^2)& 2(q_2{q}_3+q_0{q}_1)\\ 2(q_1{q}_3+q_0{q}_2)& 2(q_2{q}_3-q_0{q}_1)& 1-2(q_1^2+q_2^2)\end{array}\right)---\left(8\right)$

上式的矢量形式为$C_o^b=(q_0^2-{\overrightarrow{q}}^T\overrightarrow{q})I_{3\times 3}+2\overrightarrow{q}{\overrightarrow{q}}^T-2q_{0}S\left(\overrightarrow{q}\right).$以轨道系为参考坐标系，则由四元数来描述得到的卫星姿态运动学方程为

${\stackrel{.}{q}}_0=-1/2{\overrightarrow{q}}^T{\overrightarrow{\omega }}_{\mathrm{ob}}^b,$$\stackrel{\stackrel{.}{\to }}{q}=1/2[q_0{I}_{3\times 3}+S\left(\overrightarrow{q}\right)]{\overrightarrow{\omega }}_{\mathrm{ob}}^b---\left(9\right)$

其中，为星体系相对轨道系的角速度。相对姿态控制时，记实际四元数q与期望的四元数q_d的误差为角速度误差为则其相互关系为$\tilde{q}=q_d^{-1}\otimes q$及$\stackrel{\tilde{\to }}{\omega }={\overrightarrow{\omega }}_{\mathrm{db}}^b={\overrightarrow{\omega }}_{\mathrm{ob}}^b-R_d^b{\overrightarrow{\omega }}_{\mathrm{od}}^d={\overrightarrow{\omega }}_{\mathrm{ob}}^b-{\overrightarrow{\omega }}_{\mathrm{od}}^b,$进而得到姿态四元数误差运动方程为

${\stackrel{\stackrel{.}{~}}{q}}_0=-1/2{\stackrel{\tilde{\to }}{q}}^T\stackrel{\tilde{\to }}{\omega ,}$$\stackrel{\stackrel{\stackrel{.}{~}}{\to }}{q}=1/2[{\tilde{q}}_0{I}_{3\times 3}+S\left(\stackrel{\tilde{\to }}{q}\right)]\stackrel{\tilde{\to }}{\omega }---\left(10\right)$

在气动力矩控制下，卫星的动力学方程为

$I{\stackrel{\stackrel{.}{\to }}{\omega }}_{\mathrm{ib}}^b+{\overrightarrow{\omega }}_{\mathrm{ib}}^b\times \left(I{\overrightarrow{\omega }}_{\mathrm{ib}}^b\right)={\overrightarrow{T}}^b---\left(11\right)$

其中，I为卫星的惯量矩阵，为星体系相对于惯性系的角速度，为作用于星体的气动力矩，即${\overrightarrow{T}}^b=C_o^b{\overrightarrow{T}}_a^o.$与的关系为${\overrightarrow{\omega }}_{\mathrm{ob}}^b={\overrightarrow{\omega }}_{\mathrm{ib}}^b-C_o^b{\overrightarrow{\omega }}_{\mathrm{io}}^o={\overrightarrow{\omega }}_{\mathrm{ib}}^b+{\omega }_o{\overrightarrow{c}}_2,$其中，为式(8)中的第2列向量，ω_o为轨道角速度。

姿态确定方法主要有两个分支，其一为最小二乘估计问题，其二为估计滤波法。最小二乘法通常称为Wahba问题：在参考系中，有一组参考矢量在星体系中对应另一组观测矢量其中n≥2，根据这两组矢量来求解两坐标系之间的姿态矩阵C，使得如下误差能量函数最小

$L\left(C\right)=\underset{i=1}{\overset{n}{\Sigma }}{\sigma }_i{[{\overrightarrow{b}}_i-C{\overrightarrow{r}}_i]}^T[{\overrightarrow{b}}_i-C{\overrightarrow{r}}_i]---\left(12\right)$

其中σ_i为权系数，函数还可以写为$L\left(C\right)=\underset{i=1}{\overset{n}{\Sigma }}{\sigma }_i-\mathrm{trace}\left({\mathrm{CB}}^T\right),$$B=\underset{i=1}{\overset{n}{\Sigma }}{\sigma }_i{\overrightarrow{b}}_i{\left({\overrightarrow{r}}_i\right)}^T.$则Wahba最小能量函数问题可以转换为求如下增益函数的最大值

g(C)＝trace(CB^T)                                (13)

姿态确定的估计滤波方法采用EKF，其递推方程组如下

$P_{k/k-1}={\phi }_{k/k-1}{P}_{k-1}{\phi }_{k/k-1}^T+{\Gamma }_{k-1}{Q}_{k-1}{\Gamma }_{k-1}^T---\left(15\right)$

$K_k=P_{k/k-1}{H}_k^T{[H_k{P}_{k/k-1}{H}_k^T+R_k]}^{-1}---\left(16\right)$

$P_k=(I-K_k{H}_k)P_{k/k-1}{(I-K_k{H}_k)}^T+K_k{R}_k{K}_k^T---\left(18\right)$

其中，K_k为滤波增益，P_k/k-1和P_k分别为一步预测和估计均方差。式(14～15)为时间更新方程，式(16～18)为量测更新方程。

本发明的基于大气阻力的微小卫星编队飞行控制系统中，队形控制和姿态确定都应用到轨道信息，采用GPS进行轨道确定，当没有GPS信号时采用递推滤波进行轨道估计。当处于阴影区时姿态确定采用纯磁测加卫星动力学方程进行无陀螺估计；当处于光照区时应用太阳电池阵列和磁强计测量进行高精度融合估计。轨道确定和纯磁测姿态确定都采用EKF方法，滤波步骤如式(14～18)，下面重点对本发明所设计的双重姿态融合滤波器进行说明。

本发明的双重滤波器设计首先采用q-method法进行双矢量定姿，它是采用四元数来求解Wahba问题的有效方法。将姿态矩阵(8)代入式(13)得

$g\left(q\right)=(q_0^2-{\overrightarrow{q}}^T\overrightarrow{q})\mathrm{trace}\left(B\right)+2\mathrm{trace}\left(\overrightarrow{q}{\overrightarrow{q}}^T{B}^T\right)+2q_0\mathrm{trace}\left(\overrightarrow{q}{B}^T\right)=q^T\mathrm{Kq}---\left(19\right)$

其中，K为4×4维矩阵，定义为

$K=\left(\begin{array}{cc}B+B^T-\mathrm{trace}\left(B\right)I_{3\times 3}& Z\\ Z^T& \mathrm{trace}\left(B\right)\end{array}\right)---\left(20\right)$

其中，Z＝[B₂₃-B₃₂ B₃₁-B₁₃ B₁₂-B₂₁]^T。则K的最大特征值λ_max对应的特征向量即为所求的最优四元数q_opt，即Kq_opt＝λ_maxq_opt。

当应用q-method算法得到四元数后，进一步设计EKF滤波器。由于四元数存在约束条件‖q‖＝1，如果直接用四元数作为状态量，EKF协方差矩阵会出现奇异。本发明的滤波器设计中采用伪测量方法，在滤波的状态变量中除去四元数标量部分，且在滤波器状态更新时将四元数与其它状态量分别计算，四元数采用乘法而角速度采用加法运算，然后再统一进行滤波循环，该方法不仅可以保持四元数的正规化，而且降低状态维数从而降低计算量。观测量为由q-method计算得到四元数q_k的矢量部分，此时观测方程化为线性方程，滤波增益K_k为6×3维矩阵，由四元数和角速度两部分组成即K_k＝[K_q K_ω]^T。状态估值计算为

${\hat{q}}_k=\left(\begin{array}{c}\sqrt{1-{\left|\right|{\overrightarrow{q}}_{k/k-1}\left|\right|}^2}\\ {\overrightarrow{q}}_{k/k-1}\end{array}\right)\otimes \left(\begin{array}{c}\sqrt{1-{\left|\right|K_q{\overrightarrow{y}}_k\left|\right|}^2}\\ K_q{\overrightarrow{y}}_k\end{array}\right)---\left(21\right)$

${\stackrel{\hat{\to }}{\omega }}_{\mathrm{ob},k}^b={\overrightarrow{\omega }}_{\mathrm{ob},k/k-1}^b+K_{\omega }(q_k\otimes q_{k/k-1}^{-1})---\left(22\right)$

对于无陀螺EKF滤波器设计，上述的四元数正规化处理方法，状态量选为四元数矢量部分和角速度的6×1维矢量$X={\left(\begin{array}{cccccc}{q}_1& q_2& q_3& {\omega }_{\mathrm{ob},x}^b& {\omega }_{\mathrm{ob},y}^b& {\omega }_{\mathrm{ob},z}^b\end{array}\right)}^T,$则由卫星运动学与动力学模型，可得非线性状态方程为

$\stackrel{.}{X}=\left(\begin{array}{c}\stackrel{\stackrel{.}{\to }}{q}\\ {\stackrel{\stackrel{.}{\to }}{\omega }}_{\mathrm{ob}}^b\end{array}\right)=\left(\begin{array}{c}1/2[q_0{I}_{3\times 3}+S\left(\overrightarrow{q}\right)]{\overrightarrow{\omega }}_{\mathrm{ob}}^b\\ I^{-1}[-({\overrightarrow{\omega }}_{\mathrm{ob}}^b+C_o^b{\overrightarrow{\omega }}_{\mathrm{io}}^o)\times I({\overrightarrow{\omega }}_{\mathrm{ob}}^b+C_o^b{\overrightarrow{\omega }}_{\mathrm{io}}^o)+{\overrightarrow{T}}^b]+S\left({\overrightarrow{\omega }}_{\mathrm{ob}}^b\right)C_o^b{\overrightarrow{\omega }}_{\mathrm{io}}^o\end{array}\right)---\left(23\right)$

其中，S(·)为斜对称矩阵。对状态方程围绕最优估值进行线性化得

$F={\left(\begin{array}{c}{F}_{\overrightarrow{q}}\\ F_{\overrightarrow{\omega }}\end{array}\right)}_{X={\hat{X}}_k}=\frac{\partial f\left[X_k\right]}{\partial X_k^T}{|}_{X={\hat{X}}_k}={\left(\begin{array}{cccc}\frac{\partial \stackrel{\stackrel{.}{\to }}{q}}{\partial q_1}& \frac{\partial \stackrel{\stackrel{.}{\to }}{q}}{\partial q_2}& ···& \frac{\partial \stackrel{\stackrel{.}{\to }}{q}}{\partial {\omega }_{\mathrm{ob},z}^b}\\ \frac{\partial {\overrightarrow{\omega }}_{\mathrm{ob}}^b}{\partial q_1}& \frac{\partial {\overrightarrow{\omega }}_{\mathrm{ob}}^b}{\partial q_2}& ···& \frac{\partial {\overrightarrow{\omega }}_{\mathrm{ob}}^b}{\partial {\omega }_{\mathrm{ob},z}^b}\end{array}\right)}_{X={\hat{X}}_k}---\left(24\right)$

对于太阳电池阵列和磁强计的测量，最直接的观测量为太阳矢量和磁场矢量的量测值，即$Z={\left(\begin{array}{cccccc}{s}_{m,x}^b& s_{m,y}^b& s_{m,z}^b& b_{m,x}^b& b_{m,y}^b& b_{m,z}^b\end{array}\right)}^T,$则量测方程为

$Z=h\left[X\right]+V=\left(\begin{array}{c}{C}_o^b\left(q\right){\overrightarrow{s}}^o\\ C_o^b\left(q\right){\overrightarrow{b}}^o\end{array}\right)+V---\left(25\right)$

此时观测方程为非线性，在滤波每一步中都需要由滤波估计值得到的进行坐标转换，输出方程非常复杂。采用本发明的q-method算法进行预处理，作为EKF滤波循环外的一部分，这样可以将6维观测量转换为3维观测量，滤波循环中不需要进行坐标转换，同时观测方程也化成线性方程，大大减少了计算量。记q-method法计算得到四元数为q_c，选取q_c的矢量部分为观测量，则观测量为Z＝[q_c，x q_c，y q_c，z]^T，观测方程转化为线性方程

$Z=\mathrm{HX}+V=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)X+V_c---\left(26\right)$

其中，V_c为量测噪声，即使用q-method的计算误差。根据以上的状态方程和预处理观测方程，带入EKF递推方程组(14～18)即可进行姿态估计。

4、基于大气阻力的卫星编队飞行的姿态控制律设计：

卫星的姿态控制通常需要姿态和角速率信息，本发明可以在没有角速率测量单机(如陀螺)的前提下实现卫星的状态控制，有利于微小卫星的系统简化和节省资源。本发明的编队卫星姿态控制主要分为消除星箭分离干扰的速率阻尼和进行姿态捕获建立三轴稳定的两个阶段。

首先对于速率阻尼阶段，记期望气动板产生的控制力矩为设计控制律

${\overrightarrow{T}}_d^b=-K_{\det }{\stackrel{\stackrel{.}{\to }}{B}}^b---\left(27\right)$

其中，K_det为正定控制增益，为简便期间，选为对角阵K_det＝diag{[k_x k_y k_z]}。为磁场变化率，由磁强计的测量值进行差分计算得到的估计值代替。该方法借鉴于目前国际上使用磁控进行速率阻尼的B-dot控制方法，所不同的是控制机构为气动板而不是磁力矩器。可以近似等于卫星的角速率变化，由于磁强计所测磁场变化是由卫星飞行引起和卫星状态引起的地磁变化之和，由于飞行位置变化引起一个控制周期的磁场变化极小，因此主要反映的姿态角速度变化，根据的信息进行控制即可以阻尼角速率。

对于姿态捕获建立三轴稳定阶段，设计限制姿态反馈控制律

${\overrightarrow{T}}_d^b=-K_p·L_p\left(\overrightarrow{q}\right)-K_d{\overrightarrow{\omega }}_{\mathrm{ob}}^b---\left(28\right)$

$L_p\left(q_i\right)=\left(\begin{array}{c}{q}_i,\left|q_i\right|<L_{p\max ,i}\\ \mathrm{sign}\left(q_i\right)L_{p,\max },\left|q_i\right|\ge L_{p\max ,i}\end{array}\right),i=x,y,z$

其中，L_pmax为姿态反馈限制系数，可以通过调整L_pmax来保证该控制律满足大角度姿态捕获的要求，K_p和K_d为正定反馈增益。该控制律是在PD控制的基础上，采用限制姿态反馈的方法，使得在大角度姿态偏差时，角速度控制起主要作用，当姿态控制到一定范围后进入角速度和姿态共同作用，即为小角度下的PD控制，适合整个姿态捕获和三轴稳定阶段。

当设计好速率阻尼和姿态捕获及三轴稳定的控制律之后，需要根据期望的力矩来控制4块气动板伸缩的长度d_i(i＝1，2，3，4)，从而产生需要的气动力矩气动力矩是卫星姿态q和气动板长度d_i的非线性函数，当卫星姿态和气动板长度已知时，由式(7)可以计算控制力矩，然而由期望的控制力矩求解气动板长度的逆运算将比较困难，本发明采用线性化处理的逆运算了来设计气动板驱动控制律。气动板的控制，实际产生的气动力矩与期望的控制力矩误差为${\overrightarrow{T}}_e^b={\overrightarrow{T}}_c^b-{\overrightarrow{T}}_d^b,$在一个控制周期内期望力矩为常数，因此可得到${\stackrel{\stackrel{.}{\to }}{T}}_e^b={\stackrel{\stackrel{.}{\to }}{T}}_c^b=A{\stackrel{\stackrel{.}{\to }}{d}}_c,$其中为4块气动板的长度变化率，A为雅可比矩阵

设计气动板驱动控制律时设定初始的气动板长度为则采用气动板控制律为

其中，K_f和K_p为正定增益矩阵，均取为对角阵，为A的逆矩阵。

5、基于大气阻力的卫星编队飞行的队形控制律设计：

参见附图4本发明的基于大气阻力微小卫星编队飞行控制策略，当姿态控制完成后，通过姿态信息判断指向正确后进行编队的队形控制。编队队形控制是应用多颗卫星间的差分气动阻力进行控制：为了保持相对姿态，按同样的比例同时伸长某一卫星的4块气动板，来改变其迎风面，从而改变卫星间的气动阻力，导致卫星相对于参考坐标产生不同的加速度，进而改变卫星的相对位置。

采用Hill方程(亦称C-W方程)描述编队卫星的相对运动，根据国内文献资料的习惯重新定义轨道坐标系：x轴指向飞行方向，y轴背向地心方向，z轴垂直轨道面与x和y构成右手坐标，则Hill方程为

$\stackrel{..}{x}-2{\omega }_o\stackrel{.}{y}=f_x/m,$$\stackrel{..}{y}+2{\omega }_o\stackrel{.}{x}-3{\omega }_o^{2}y=f_y/m,$$\stackrel{..}{z}+{\omega }_o^{2}z=f_z/m---\left(31\right)$

其中，ω_o为轨道角速度，$\overrightarrow{F}={\left(\begin{array}{ccc}{f}_x& f_y& f_z\end{array}\right)}^T$为相对于卫星各轴的作用力，m为卫星的质量。当没有作用力即编队卫星自由飞行时，对Hill方程进行求解分析，可以将卫星运动分解为轨道面(xy平面)和垂直于轨道面(z方向)的两个独立运动，xy平面内相对运动为椭圆，z方向为简谐振动。

大气阻尼会造成卫星编队构形产生沿航迹方向的漂移，而不产生垂直于轨道面的作用力，因此本发明的编队队形控制仅考虑共面的编队飞行。卫星的编队构形主要有两类：其一为主从结构编队(leader-follower)，其二为虚拟参考卫星编队。主从结构编队构型的主卫星位于参考坐标系，从卫星相对于主卫星运动；虚拟参考卫星编队的参考中心没有实质的卫星，各卫星相对于某个参考中心运动，是一种更为稳定的编队形式，不会因为主卫星失效而使整个编队系统实效。因此，下面主要针对共面虚拟参考卫星编队的两种编队队形进行控制律设计：共面跟飞和共面椭圆编队，参见附图4a、附图4b的基于大气阻力的微小卫星编队飞行队形示意图。

对于共面编队卫星，以双星编队为例设计控制律，选取控制状态量为

$X\left(t\right)={\left(\begin{array}{cccccccc}{x}_1& y_1& {\stackrel{.}{x}}_1& {\stackrel{.}{y}}_1& x_2& y_2& {\stackrel{.}{x}}_2& {\stackrel{.}{y}}_2\end{array}\right)}^T---\left(32\right)$

大气阻力在卫星各轴向的分量为$f_x={-c}_d\rho V_R^{2}A/2,$f_y＝0，$f_z={\mathrm{ac}}_d\rho V_R^{2}A/2,$由于a＝(1.5ω_e/ω_o)sin i cos ω_ot<1×10^-4，则f_z小于f_x的10^-4量级，故可不做考虑。记c_dA_i/m_i＝γ_i(i＝1，2)，选取控制量为u＝[γ₁ γ₂]^T，则系统状态方程为

$\stackrel{.}{X}\left(t\right)=\mathrm{AX}\left(t\right)+\mathrm{Bu}\left(t\right)---\left(33\right)$

$A={\left(\begin{array}{cc}{A}_1& 0_{4\times 4}\\ 0_{4\times 4}& A_2\end{array}\right)}_{8\times 8},$

$B={\left(\begin{array}{cc}{B}_1& 0_{4\times 1}\\ 0_{4\times 1}& B_2\end{array}\right)}_{8\times 2}$

$A_1=A_2=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 3{\omega }_o^2& 0& 0& -2{\omega }_o\\ 0& 0& 2{\omega }_o& 0\end{array}\right),$$B_1=B_2={\left(\begin{array}{cccc}0& 0& 0& -\rho V_R^2/2\end{array}\right)}^T$

本发明的卫星编队队形控制采用最优控制方法的LQ终端控制器进行控制律设计。采用大气阻力控制卫星时，主要通过伸缩气动板长度来改变迎风面面积，从而改变控制量u＝[γ₁ γ₂]^T。但是大气阻力增加将降低卫星速度，进而降低卫星轨道半径，因此编队控制的同时应尽量减少轨道降低。定义轨道能量函数为$ϵ=V_R^2/2-\mu /r=-\mu /2a,$其中r为轨道速度和半径，a为轨道半长轴，则减少轨道半径将降低轨道能量。对于卫星编队系统(33)，设系统初始状态为X(t₀)＝X₀，采用LQ软终端约束控制的目标是极小化指标泛函

$J=\frac{1}{2}{\int }_0^{t_f}\left(u^T\mathrm{Ru}\right)\mathrm{dt}+\frac{1}{2}{e}_f^T{Q}_f{e}_f=\frac{1}{2}{\int }_0^{t_f}({\left|\frac{{\mathrm{dϵ}}_1}{\mathrm{dt}}\right|}^2+{\left|\frac{{\mathrm{dϵ}}_2}{\mathrm{dt}}\right|}^2)\mathrm{dt}+\frac{1}{2}{(M_{f}X\left(t_f\right)-\psi )}^T{Q}_f(M_{f}X\left(t_f\right)-\psi )$

(34)

其中，u＝[c_dA₁/m₁ c_dA₂/m₂]^T为控制量，$R={\left(\begin{array}{ccc}1/4{\rho }^2{V}_R^6& 0;0& 1/4{\rho }^2{V}_R^2\end{array}\right)}^T$为正定实对称矩阵；Q_f为半正定实对称矩阵，e_f为终端误差向量，M_f和ψ分别为终端状态矩阵和向量。由系统状态方程和指标泛函，得哈密顿函数

H＝1/2u^T Ru+λ^T(AX+Bu)                   (35)

其中，λ为Lagrange乘子，则由最优性的必要条件得到最优控制律u(t)为

u(t)＝-R^-1B^Tλ(t)                        (36)

由控制律(36)及哈密顿函数(35)可以得到哈密顿正则微分方程，以及软终端约束LQ控制导出的哈密顿两端边值问题的边界条件分别为。

$\left(\begin{array}{c}\stackrel{.}{X}\left(t\right)\\ \stackrel{.}{\lambda }\left(t\right)\end{array}\right)=\left(\begin{array}{cc}A& -{\mathrm{BR}}^{-1}{B}^T\\ 0_{8\times 8}& -A^T\end{array}\right)\left(\begin{array}{c}X\left(t\right)\\ \lambda \left(t\right)\end{array}\right)---\left(37\right)$

X(t₀)＝X₀，$\lambda \left(t_f\right)=M_f^T{Q}_f[M_{f}X\left(t_f\right)-\psi ]---\left(38\right)$

不同的卫星编队队形具有不同的终端约束，通过对本发明附图5的两种编队卫星任务及队形特征的分析，得到共面跟飞的终端状态矩阵和向量为

$M_f=\left(\begin{array}{cc}{M}_1& 0_{4\times 4}\\ 0_{4\times 4}& M_2\end{array}\right),$$M_1=\left(\begin{array}{cccc}2{\omega }_0& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right),$$M_1=\left(\begin{array}{cccc}-2{\omega }_0& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\end{array}\right),$

ψ＝[0 l/2 0 0 0 l/2 0 0]^T                                       (39)

共面椭圆编队的终端矩阵M_f与共面跟飞编队的相同，所不同的是终端向量，为

ψ＝[0 l ω_ol/2 0 0 l ω_ol/2 0]^T                             (40)

对于本发明的基于大气阻力的共面编队飞行队形控制，根据编队任务的要求，如队形建立、队形重构以及队形保持等任务，设置初始条件和终端目标后，便可以根据上述的LQ终端控制律计算得到控制量u＝[γ₁ γ₂]^T，即可以得到控制量为[c_dA₁/m₁ c_dA₂/m₂]^T。卫星的阻力系数c_d和质量m_i(i＝1，2)是固定的，则由控制量u计算需要的迎风面A，而迎风面大小可以由气动板的伸缩来控制。

6、有益效果分析：

根据以上五个主要技术实施方式，对本发明的有益效果进行分析。首先针对步骤2所设计的编队卫星建立卫星相对运动模型，在此基础上对编队飞行姿态确定与控制、编队队形控制进行分析。根据目前国内外的微小卫星设计指标分析，设编队飞行微小卫星质量均为m＝20kg，转动惯量I＝[1.0 1.05 0.9]^Tkg·m²，卫星尺寸为a×b×c＝0.6×0.6×0.7m³，4块气动板的宽均为e＝0.3m，长度的控制范围为0m≤d≤0.35m，编队飞行高度为600km。

对于编队卫星飞行姿态确定，根据磁强计及太阳电池阵列测量性能指标的实际特性分析，取磁强计测量误差为300nT(3σ)，太阳电池阵列测量误差2°(3σ)，滤波周期为1s。姿态捕获阶段的姿态运动变化较为丰富，因此选取该阶段进行测量，运动轨迹如附图5a，三轴初始姿态为[40 60 150]°，采用6维直接观测量的结果、单独q-method方法测量结果和本发明的双重滤波结果如图5b、5c、5d所示，可以看出采用直接EKF滤波器时，在0～2000s大姿态变化时误差达到20°，而使用本发明的双重滤波后误差在5°以内，在2000s以后的稳态阶段，不同滤波器的比较结果见表1所示，可以看出双重滤波器具有最高的精度，而且本发明的双重滤波器将6维观测量转化为四元数，将观测方程化为线性，计算量要比直接EKF滤波器大大减小。

表1  滤波稳定后不同无陀螺滤波器结果比较

对于基于大气阻力的编队卫星姿态控制，采用前面的卫星设计指标，气动板控制变化范围为0m≤d≤0.35m，给出速率阻尼、三轴姿态控制情况的有益效果。速率阻尼时卫星三轴初始角速度为[0.6 0.6 0.6]°/s，约为10倍的轨道角速度，控制结果如附图5e、5f所示。卫星三轴姿态控制时初始姿态为[35 30 -20]°，控制目标为三轴姿态均为0°，即姿态四元数为[1 0 0 0]，控制结果如附图5g、5h所示。可以看出本发明的基于大气阻力的编队卫星姿态控制具有效果良好，通过多次实验结果表明在编队卫星处于1000km以下中、低轨道时均可有效的进行姿态控制，轨道越低则达到目标姿态的控制速度越快。

对于基于大气阻力的编队卫星队形控制，为了保持编队卫星的相对姿态，各卫星4块气动板应以相同的比例进行伸缩控制，若队形控制开始时各颗卫星的四块气动板长度相同，则队形控制过程中应保持d₁＝d₂＝d₃＝d₄。对于附图4a的共面跟飞和附图4b共面椭圆两种队形控制，采用本发明所设计的队形控制律进行控制，队形建立时设两星在x和y方向的初始位置为[-10 10]m和[-10 10]m，初始速度为[-0.02 0.02]m/s和[0.02 -0.02]m/s，共面跟飞终端目标为两星位于同一轨道相距1000m，椭圆编队终端目标为半长轴500m，控制时间为5个轨道周期，结果如附图5i、5j，其中符号“”为编队相对参考中心；队形重构控制时，共面跟飞→椭圆编队重构的初始距离为1000m，终端目标为椭圆半长轴100m，共面椭圆→跟飞编队重构时初始椭圆长半轴500m，终端目标相距100m，结果如附图5k、5l所示。可以看出本发明的基于大气阻力的编队卫星队形控制具有效果良好，通过多次实验结果表明在编队卫星处于1000km以下中、低轨道时均可有效的进行队形建立、重构以及保持等控制。