立方星陀螺/双星敏感器组合定姿方法*

马海宁，陆正亮，张 翔，廖文和

(南京理工大学 机械工程学院， 江苏 南京 210094)

立方星具有质量轻、体积小、制造周期短、功能密度高等特点，在军用和民用领域具有广阔的应用前景[1-2]。姿态确定系统提供实时的姿态信息，其精度直接影响立方星的姿态控制精度，从而影响卫星的工作性能[3]。陀螺/星敏感器组合的姿态确定系统是现代高精度卫星的主要定姿方式，随着低成本、小型化的微机电系统(Micro-Electro-Mechanical System, MEMS)陀螺和纳型星敏感器的不断发展，其在立方星上的应用正逐步普及[4-5]。

近年来，国内外很多学者基于陀螺/星敏感器的配置，采用不同的定姿方法对卫星的姿态确定进行了深入研究。文献[6-8]基于陀螺和星敏感器，利用扩展卡尔曼滤波(Extended Kalman Filter, EKF)算法估计卫星的姿态并对陀螺偏置进行补偿。文献[9]基于陀螺和星敏感器，采用实时的滚动时域估计方法获取卫星姿态和陀螺校准参数，该方法对于较大的初始姿态误差具有更快的收敛速度。文献[10]基于陀螺和星敏感器，提出利用容积卡尔曼滤波(Cubature Kalman Filter, CKF)算法估计卫星姿态，结果表明该方法在处理姿态确定系统的强非线性问题时比EKF算法和无迹卡尔曼滤波(Unscented Kalman Filter, UKF)算法拥有更高的精度和稳定性。文献[11]针对陀螺和星敏感器噪声频率特性不同的特点，提出利用互补滤波算法获得姿态融合数据，将其引入UKF算法的测量更新过程，最终获得卫星姿态的最优估计。

上述研究均是针对单星敏感器配置的情况，然而由于星敏感器视场角较小，导致姿态解算时滚转轴的精度低于俯仰轴和偏航轴精度1～2个数量级[12]。因此，有关学者开展了多星敏感器姿态测量的研究：文献[13]提出利用基于UKF的联邦滤波算法融合双星敏感器的测量数据以实现卫星的姿态估计，并采用故障检测和隔离算法规避敏感器测量故障，但实际上仍是主要利用主星敏感器的测量数据，副星敏感器只是作为备份；文献[14]提出了基于星敏感器观测姿态角的误差建模，研究了多星敏感器组合的最优安装构型，利用陀螺加星敏感器组合的EKF算法对卫星姿态和陀螺漂移进行估计，但姿态估计精度不高；文献[15]提出利用基于EKF的联邦滤波算法融合三星敏感器的测量数据以实现卫星的姿态估计，该方法提高了定姿精度和系统鲁棒性，保证了定姿任务的持续性，但其采用的星敏感器是应用于传统大卫星的高精度、昂贵的大型星敏感器，不适用于低成本、低功耗要求的立方星平台。

鉴于此，本文结合某型号立方星研制任务的需要，提出了陀螺/双星敏感器的组合定姿方法。该方法在硬件层面，采用两个同时工作的低成本、低精度的皮型星敏感器以及三轴MEMS陀螺作为姿态敏感器；在算法层面，基于平均四元数的思想，构建了两种姿态估计方案：一种是先对双星敏感器输出的姿态四元数进行平均，再结合陀螺的测量，利用传统的MEKF算法实现卫星姿态和陀螺漂移的最优估计，即集中式定姿；另一种是先利用传统的MEKF算法分别融合陀螺/星敏感器1和陀螺/星敏感器2的测量，再对两者获得的姿态四元数和陀螺漂移进行平均，得到卫星姿态和陀螺漂移的最优估计，即分散式定姿。最后通过数值仿真实验对所提出的定姿方法进行了验证。

1 姿态运动学和敏感器模型

1.1 姿态运动学模型

由于四元数具有非奇异性，因此采用四元数来描述卫星的姿态，定义姿态四元数为：

(1)

用姿态四元数表示的卫星姿态运动学方程为：

(2)

式中，qbo为轨道系相对于卫星本体系的姿态四元数，ωbo为卫星本体系相对于轨道系的旋转角速度在卫星本体系中的表达，Ω(ωbo)为四元数右乘矩阵，即

(3)

1.2 MEMS陀螺测量模型

MEMS陀螺三个测量轴分别对准卫星本体系的惯量主轴三个方向安装，它能提供卫星本体系相对于惯性系的角速度信息。考虑陀螺漂移、随机误差和测量噪声等因素的存在，构建离散测量模型为：

(4)

1.3 星敏感器

1.3.1 安装方位分析

本研究选用的星敏感器为国内某公司生产的型号为PST-3的皮型星敏感器，其杂光抑制角θl<35°。对于双星敏感器配置，为了获得较高的测量精度以及完整的可观测性，通常呈垂直关系安装[16]。立方星本体共有6个体装面，分别为对天面、对地面、迎风面、背风面、向阳面和背阴面。星敏感器在轨运行时主要受到太阳光和地气光的干扰，其中太阳光主要照射面为对天面、对地面、迎风面、背风面和向阳面，次要照射面为背阴面；地气光照射面为对地面。因此，考虑将双星敏感器垂直安装于卫星本体的背阴面，并通过系统工具包(Systems Tool Kits, STK)的建模仿真验证此种安装方式能够有效避免太阳光对星敏感器视场的干扰。图1为在STK中建立的卫星在轨飞行模型。

图1 卫星在轨飞行模型Fig.1 Flight model of satellite in orbit

卫星进入太阳阴影区后，星敏感器不受太阳光的照射，故视场分析不必考虑该区间。根据卫星的轨道参数和初始历元时间2018年1月19日12:13:39.790 UTCG，利用STK获取其在轨运行一年内地日连线与轨道平面的夹角(即太阳入射角)的变化情况，如图2所示。

图2 卫星在轨一年内太阳入射角Fig.2 Solar incident angles of satellite within one year in orbit

由图2可见，太阳入射角在第3192 h达到最小，对应的日期约为2018年6月1日12时15分。此时星敏感器视轴和太阳矢量的夹角最小，太阳光对星敏感器的影响最大。因此，本文在分析太阳光对星敏感器视场的影响时仅考虑光照区最为恶劣的情况，仿真时间设置为2018年6月1日12:13:52.862 UTCG至2018年6月1日13:14:17.117 UTCG。图3所示为上述光照区内，星敏感器视轴与太阳矢量夹角的变化情况。

图3 星敏感器视轴与太阳矢量夹角Fig.3 Angles between star sensors and sun vector

由图3可见，2个星敏感器与太阳矢量的最小夹角均大于杂光抑制角θl，这说明本文双星敏感器的安装方式使得卫星整个在轨寿命期间能有效地避免太阳光对星敏感器视场的干扰。

1.3.2 测量模型

双星敏感器的安装示意如图4所示。

图4 双星敏感器安装示意图Fig.4 Double star sensors installation diagram

单星敏感器能测量出惯性系相对于星敏本体系的姿态四元数，其测量模型为：

(5)

(6)

(7)

式中，ψbs、φbs、θbs为星敏感器在卫星本体中的安装角度。

2 非线性滤波算法

考虑到稳定状态下立方星姿态运动学和动力学方程非线性较弱，且星载计算机运算能力有限，因此选用迭代速度快且估计精度较高的MEKF算法作为定姿滤波器。MEKF的基本思想是估计无约束的三分量姿态误差参数并利用四元数乘法为航天器提供全局非奇异姿态描述[17]。

2.1 滤波状态方程

(8)

(9)

对式(9)作线性化处理，得线性连续滤波状态方程为：

(10)

对式(10)作离散化处理，得线性离散滤波状态方程为：

(11)

2.2 滤波观测方程

(12)

取出式(12)两端的矢量部分：

(13)

由式(13)构建线性离散滤波观测方程：

δZk=HkΔXk/(k-1)+Vk

(14)

2.3 姿态估计流程

2.3.1 状态量预估

k-1时刻，利用姿态运动学方程和陀螺测量值，计算姿态四元数的预测值和陀螺漂移的预测值：

(15)

(16)

2.3.2 状态偏差量更新

一步预测均方误差方程：

(17)

滤波增益方程：

(18)

状态估值方程：

(19)

最优估计均方误差方程：

(20)

2.3.3 状态量修正

(21)

3 陀螺/双星敏感器定姿

针对双星敏感器的配置，需要考虑对两者的输出四元数进行融合，以获得最终的测量数据并传入MEKF滤波器的观测方程中。由于四元数的非加性特点，常见的线性加权方法由于无法保证四元数归一化约束以及非唯一性特性，并不能直接处理多个四元数的平均。针对此问题，文献[18]给出了一种加权四元数方法；文献[19]给出了一种优良的四元数平均方法，可以直接用于平均两个四元数：

(22)

基于文献[19]的研究结果以及非线性加权平均的思想，融合双星敏感器的输出四元数，同时可延伸到对多个姿态四元数估计值的融合。由此提出两种陀螺/双星敏感器组合定姿方法，即集中式和分散式。

1)集中式：如图5所示，先利用式(22)对双星敏感器的姿态四元数测量值进行平均处理，再结合MEMS陀螺的测量，利用MEKF算法实现卫星姿态和陀螺漂移的最优估计。

图5 集中式定姿方法Fig.5 Centralized attitude determination method

2)分散式：如图6所示，先利用MEKF算法分别融合MEMS陀螺和星敏感器1以及MEMS陀螺和星敏感器2的测量数据，获得各自的姿态四元数估计和陀螺偏置估计，再利用式(22)进行四元数加权平均处理，获得姿态最优估计，同时利用线性加权平均方法获得陀螺偏置的最优估计。

图6 分散式定姿方法Fig.6 Decentralized attitude determination method

4 数值仿真及结果分析

利用数值仿真对本文提出的陀螺/双星敏感器组合定姿方法进行验证。

4.1 仿真参数设定

1) 轨道高度为500 km，轨道倾角为97.5°，轨道偏心率为0.001，轨道角速度ωo为0.001 1 rad/s。

3) 星敏感器视轴方向测量误差Δva为50″，视平面方向测量误差Δfp为5″；星敏感器1安装角度为ψbs1=-90°,φbs1=45°,θbs1=-90°，星敏感器2安装角度ψbs2=90°,φbs2=45°,θbs2=90°；星敏感器1测量权重w1=0.5，星敏感器2测量权重w2=0.5，采样频率为10 Hz。

4) 三轴姿态角初始真值为[2° -2° 2°]，三轴角速度初始真值为[-0.5 (°)/s0.5 (°)/s0.5 (°)/s]。

5) 三轴姿态角初始估计值为[0° 0° 0°]°，三轴角速度初始估计值为[0 (°)/s0 (°)/s0 (°)/s]，陀螺偏置初始估计值：[0 (°)/h0 (°)/h0 (°)/h]。

6) 数值仿真步长为1 s，仿真时长为5000 s。

7) 仿真计算机配置：Intel Core i5-4590 3.3 GHz，四核；RAM 8 GB；MATLAB R2015a。

4.2 仿真结果及分析

根据上述参数设定，分别对陀螺/星敏感器1、集中式陀螺/双星敏感器和分散式陀螺/双星敏感器组合定姿方式进行数值仿真分析。

4.2.1 定姿精度分析

陀螺/星敏感器1组合定姿得到的稳定阶段姿态参数误差变化曲线如图7所示，三个方向的姿态角误差的变化范围分别为±0.020 4°(3σ)、±0.029 7°(3σ)、±0.020 4°(3σ)；三个方向的陀螺偏置误差的变化范围分别为±0.008 1 (°)/s(3σ)、±0.007 5 (°)/s(3σ)、±0.007 5 (°)/s(3σ)。

(a) 偏航角误差(a) Error of yaw angle

(b) 滚动角误差(b) Error of roll angle

(c) 俯仰角误差(c) Error of pitch angle

(d) 陀螺X轴偏置误差(d) Error of gyro bias in X axis

(e) 陀螺Y轴偏置误差(e) Error of gyro bias in Y axis

(f) 陀螺Z轴偏置误差(f) Error of gyro bias in Z axis图7 陀螺/星敏感器1定姿结果(稳定阶段)Fig.7 Attitude determination results based on MEMS gyro and star sensor 1(stable stage)

集中式陀螺/双星敏感器组合定姿得到的稳定阶段姿态参数误差变化曲线如图8所示，三个方向的姿态角误差的变化范围分别为±0.009 0°(3σ)、±0.018 3°(3σ)、±0.018 0°(3σ)；三个方向的陀螺偏置误差的变化范围分别为±0.007 5 (°)/s(3σ)、±0.006 9 (°)/s(3σ)、±0.007 2 (°)/s(3σ)。

分散式陀螺/双星敏感器组合定姿得到的稳定阶段姿态参数误差变化曲线如图9所示，三个方向的姿态角误差的变化范围分别为±0.012 0°(3σ)、±0.018 9°(3σ)、±0.018 6°(3σ)；三个方向的陀螺偏置误差的变化范围分别为±0.007 5 (°)/s(3σ)、±0.007 2 (°)/s(3σ)、±0.007 2 (°)/s(3σ)。

(a) 偏航角误差(a) Error of yaw angle

(b) 滚动角误差(b) Error of roll angle

(c) 俯仰角误差(c) Error of pitch angle

(d) 陀螺X轴偏置误差(d) Error of gyro bias in X axis

(e) 陀螺Y轴偏置误差(e) Error of gyro bias in Y axis

(f) 陀螺Z轴偏置误差(f) Error of gyro bias in Z axis图8 陀螺/双星敏感器集中式滤波定姿结果(稳定阶段)Fig.8 Attitude determination results of centralized filter based on MEMS gyro and double star sensors(stable stage)

表1列出了三种组合定姿方式得到的姿态角均方根误差(Root Mean Square Error, RMSE)和平均绝对误差(Mean Absolute Error, MAE)的统计结果。由表1可见，在定姿精度方面，陀螺/双星敏感器集中式定姿得到的三个方向的姿态角RMSE值和MAE值均比分散式定姿得到的对应结果小，这表明集中式定姿的精度高于分散式定姿；集中式定姿和分散式定姿得到的三个方向的姿态角RMSE值和MAE值均比陀螺/星敏感器1定姿得到的对应结果小，这表明双星敏感器定姿精度优于单星敏感器。

(a) 偏航角误差(a) Error of yaw angle

(b) 滚动角误差(b) Error of roll angle

(c) 俯仰角误差(c) Error of pitch angle

(d) 陀螺X轴偏置误差(d) Error of gyro bias in X axis

(e) 陀螺Y轴偏置误差(e) Error of gyro bias in Y axis

(f) 陀螺Z轴偏置误差(f) Error of gyro bias in Z axis图9 陀螺/双星敏感器分散式滤波定姿结果(稳定阶段)Fig.9 Attitude determination results of decentralized filter based on MEMS gyro and double star sensors(stable stage)

表1 姿态角RMSE和MAE的统计结果比较

表2列出了三种组合定姿方式得到的陀螺偏置RMSE和MAE的统计结果。由表2可见，在陀螺偏置估计精度方面，陀螺/双星敏感器集中式定姿得到的三个方向的陀螺偏置RMSE值和MAE值与分散式定姿得到的对应结果基本一致，这表明两者对陀螺偏置的估计精度相当；集中式定姿和分散式定姿得到的三个方向的陀螺偏置RMSE值和MAE值均比陀螺/星敏感器1定姿得到的对应结果略小，这表明双星敏感器对陀螺偏置的估计精度略优于单星敏感器。

表2 陀螺偏置RMSE和MAE的统计结果比较Tab.2 Statistical results comparison of the RMSE and MAE of gyro bias

4.2.2 滤波收敛性分析

陀螺/星敏感器1组合定姿得到的初始阶段姿态角误差变化曲线如图10所示，从图中可以看出，滤波收敛时间约为23 s。

集中式陀螺/双星敏感器组合定姿得到的初始阶段姿态角误差变化曲线如图11所示，从图中可以看出，滤波收敛时间约为20 s。

(a) 偏航角误差(a) Error of yaw angle

(b) 滚动角误差(b) Error of roll angle

(c) 俯仰角误差(c) Error of pitch angle图10 陀螺/星敏感器1定姿结果(初始阶段)Fig.10 Attitude determination results based on MEMS gyro and star sensor 1(initial stage)

(a) 偏航角误差(a) Error of yaw angle

(b) 滚动角误差(b) Error of roll angle

(c) 俯仰角误差(c) Error of pitch angle图11 陀螺/双星敏感器集中式滤波定姿结果(初始阶段)Fig.11 Attitude determination results of centralized filter based on MEMS gyro and double star sensors(initial stage)

分散式陀螺/双星敏感器组合定姿得到的稳定阶段姿态角误差变化曲线如图12所示，从图中可以看出，滤波收敛时间约为15 s。

由以上仿真数据可以看出，陀螺/双星敏感器集中式定姿的收敛速度与陀螺/星敏感器1定姿相当，而分散式定姿的收敛速度优于二者。

(a) 偏航角误差(a) Error of yaw angle

(b) 滚动角误差(b) Error of roll angle

(c) 俯仰角误差(c) Error of pitch angle图12 陀螺/双星敏感器分散式滤波定姿结果(初始阶段)Fig.12 Attitude determination results of decentralized filter based on MEMS gyro and double star sensors(initial stage)

5 结论

本文提出了一种陀螺/双星敏感器组合的立方星姿态确定方法，采用低成本、低精度器件作为姿态敏感器，并基于平均四元数的思想，构建了两种定姿方案：一是先对双星敏感器输出的姿态四元数进行平均，再结合陀螺的测量，利用MEKF算法实现姿态最优估计；二是先利用MEKF算法分别获取陀螺/星敏感器1和陀螺/星敏感器2的姿态估计，再对得到的姿态四元数进行平均，实现姿态最优估计。数学仿真结果表明，两种定姿方案能够提高立方星的定姿性能，其中集中式定姿获得的姿态确定精度最高，而分散式定姿的收敛速度最快。因此，本文提出的定姿方法能够为立方星实现低成本、高精度、高可靠性的姿态确定提供一种可行的工程应用参考方案。