CN110285815A - 一种可在轨全程应用的微纳卫星多源信息姿态确定方法 - Google Patents

Abstract

本发明提供一种可在轨全程应用的微纳卫星多源信息姿态确定方法，包括如下步骤：计算信息分配因子α_i，并同时向局部滤波器S1和局部滤波器S2自适应重置过程噪声方差矩阵Q_i和滤波过程估计误差协方差矩阵P_k‑1,i，下角标i＝1,2分别对应两个局部滤波器；所述信息分配因子其中，表示k‑1时刻估计误差协方差矩阵P_i,k‑1的第n行、第m列的元素；利用自适应重置后的Q_i和P_k‑1,i实现微纳卫星多源信息姿态估计。本发明加入自适应思想构成自适应联邦卡尔曼滤波算法，实现了采用多个低成本、体积小、低功耗微型敏感器的微纳卫星的定姿需求。

Description

一种可在轨全程应用的微纳卫星多源信息姿态确定方法

技术领域

本发明涉及一种卫星姿态确定方法，具体涉及一种可在轨全程应用的微纳卫星多源信息姿态确定方法。

背景技术

姿态确定与控制系统(Attitude Determination and Control System，ADCS)是微纳卫星平台在轨正常工作的基本保障。

微纳卫星因质量、体积和功耗的限制，多采用微型化的敏感器，如MEMS陀螺、磁强计、微型太阳敏感器等。微纳卫星的姿态信息可采用惯性陀螺测量，通过积分求解卫星的旋转姿态；但在实际工程应用中，MEMS陀螺的输出信息中存在漂移且易受到环境因素干扰，其转动角速度的噪声与漂移偏差值会随着积分过程持续累积与放大，造成姿态确定精度不理想。磁强计通过测量地磁场矢量确定卫星姿态，但易受到卫星自身剩磁等因素干扰，姿态确定精度不高。太阳敏感器可用于测量太阳矢量方向，但当卫星处于阴影区时不能正常工作，也就无法利用太阳矢量确定卫星姿态。因此，针对单一敏感器无法定姿或获取卫星姿态确定精度不高的问题，为保证微纳卫星姿态确定系统的定姿精度，基于多源信息融合技术，将多个敏感器数据进行融合获取卫星的姿态信息。

多源信息融合技术多采用集中式和分布式融合结构。相比于集中式融合结构，分布式融合在容错性和运算效率方面具有显著优势。多源信息融合方法最常用的是扩展卡尔曼滤波(Extented Kalman Filter,EKF)，它是将非线性系统线性化，然后再进行卡尔曼滤波。从扩展卡尔曼滤波原理可知，该方法对非线性较弱的系统具有较好的滤波效果；但当面对强非线性系统时，该方法会因多个极值点的存在造成系统滤波发散。Carlson在卡尔曼滤波基础上提出了一种分布式融合联邦滤波算法，其各子滤波器并行运算，进行状态估计，主滤波器实现信息融合输出最终结果；融合的思想主要源于卡尔曼滤波，因此也称为分布式联邦卡尔曼滤波。联邦卡尔曼滤波具有分布式融合的优点，即较好的容错性能和较少的计算量。经查阅资料发现，倪枫等人基于多传感器设计了一种基于多传感器的微小卫星姿态确定方法，该方法中未对子滤波器的过程噪声方差矩阵进行自适应分配，而且当误差协方差矩阵初值选为零矩阵时，存在奇异问题。

发明内容

有鉴于此，本发明克服现有技术的不足之处，提供一种可在轨全程应用的微纳卫星多源信息姿态确定方法，对已有技术中的联邦滤波算法进行改进，并加入自适应思想构成自适应联邦卡尔曼滤波算法，实现了采用多个低成本、体积小、低功耗微型敏感器的微纳卫星的定姿需求，并且具有容错性好、精度高、计算量小且可在轨全程应用等优点。

本发明的技术解决方案是：

一种可在轨全程应用的微纳卫星多源信息姿态确定方法，包括如下步骤：

计算信息分配因子α_i，并同时向局部滤波器S1和局部滤波器S2自适应重置过程噪声方差矩阵Q_i和滤波过程估计误差协方差矩阵P_k-1,i，下角标i＝1,2分别对应两个局部滤波器；

所述信息分配因子

其中，表示k-1时刻估计误差协方差矩阵P_i,k-1的第n行、第m列的元素；

利用自适应重置后的Q_i和P_k-1,i实现微纳卫星多源信息姿态估计。

进一步地，本发明所述自适应重置过程噪声方差矩阵Q_i和滤波过程估计误差协方差矩阵P_k-1,i为：

P_k-1,i＝α_i(P_k/k-1+I)

Q_i＝α_iQ

其中，P_k/k-1表示k时刻主滤波器中一步预测估计误差协方差矩阵，Q为主滤波器的过程噪声方差矩阵。

进一步地，本发明所述局部滤波器S1的观测信息为：

Z_1k＝(B_b(k)-B_b(k/k-1))/2

B_b(k/k-1)＝T_bo(q_bo(k/k-1))B_o(k)

其中，B_b(k)为k时刻卫星本体坐标系下地磁场矢量测量值；B_o(k)为k时刻卫星轨道坐标系下地磁场矢量；q_bo(k/k-1)为k-1时刻到k时刻的姿态四元数的预测值；T_bo(·)为卫星轨道系到本体系的方向余弦阵。

进一步地，本发明所述局部滤波器S2的观测信息为：

Z_2k＝(S_b(k)-S_b(k/k-1))/2

S_b(k/k-1)＝T_bo(q_bo(k/k-1))S_o(k)

其中，S_b(k)为k时刻卫星本体坐标系下太阳矢量测量值；S_o(k)为k时刻卫星轨道坐标系下太阳矢量。

有益效果：

第一，本发明基于多个低成本、小体积和低功耗的微型敏感器，结合多源信息分布式融合结构实现微纳卫星的姿态确定。所提出的改进型分布式自适应联邦卡尔曼滤波算法相比集中式卡尔曼定姿方法，容错性好、系统生命力强、计算量小、工程易于实现，有助于形成一套可工程化且通用化的微纳卫星姿态确定系统。

第二，本发明的改进型分布式自适应联邦卡尔曼滤波算法针对信息分配因子的修改，避免了当误差协方差矩阵初值选为零矩阵时，存在的奇异问题，保证了滤波的快速收敛及定姿精度。

第三，本发明考虑卫星在轨时姿态确定系统中过程噪声方差矩阵随环境变化的情况，通过自适应分配因子调整局部滤波器的过程噪声方差矩阵，使得所设计的姿态确定系统更贴近在轨的真实情况，因此，所发明的微纳卫星姿态确定方案具有较高的估计精度和可靠性。

第四，本发明的改进型分布式自适应联邦卡尔曼滤波算法可应用于卫星在轨的全模式，如卫星的速率阻尼模式、太阳捕获模式、对日巡航模式及对地三轴稳定模式。

附图说明

图1为基于多源信息分布式融合的基本结构图；

图2为基于多源信息融合的姿态确定方法基本原理图；

图3为改进型分布式自适应联邦卡尔曼滤波基本原理图；

图4为微纳卫星对地三轴稳定转速率阻尼模式时，基于本发明改进型分布式自适应联邦卡尔曼滤波算法所确定的卫星的姿态信息。

图5为微纳卫星对地三轴稳定转速率阻尼模式时，基于本发明改进型分布式自适应联邦卡尔曼滤波算法所确定的卫星角速度信息。

图6为微纳卫星对地三轴稳定转太阳捕获巡航时，基于本发明改进型分布式自适应联邦卡尔曼滤波算法所确定的卫星姿态信息。

图7为微纳卫星对地三轴稳定转太阳捕获巡航时，基于本发明改进型分布式自适应联邦卡尔曼滤波算法所确定的卫星角速度信息。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述。

图1给出了基于多源信息分布式融合的基本结构图，所谓多源信息融合即是多个传感器的数据融合，因此图中数据来源是敏感器1、敏感器2、……、敏感器n。然后各个敏感器数据进行量测、滤波更新，最后进入融合中心进行数据融合。本发明基于图1的多源信息分布式融合结构，选择MEMS陀螺、磁强计和太阳敏感器作为姿态敏感器，进行多源信息融合微纳卫星姿态确定方法的设计，基本原理图如图2所示。

如图3本发明中改进型分布式自适应联邦卡尔曼滤波基本原理图可知，本发明中的改进型分布式自适应联邦滤波算法由两个局部滤波器和一个主滤波器组成，是一种融合—自适应分配—反馈结构。局部滤波器S1和局部滤波器S2分别根据反馈的姿态和自适应分配的协方差和过程噪声方差完成磁强计和太阳敏感器的估计误差协方差的预测及量测更新；主滤波器实现姿态四元数的一步预测，协方差的一步预测，根据自适应分配因子进行信息分配并反馈给两个局部滤波器，还有根据两个局部滤波器的输出信息进行信息融合，然后经过姿态及陀螺漂移修正后将最终姿态信息输出。其中两个局部滤波器的状态变量相同，二者进行并行运算。

其中，计算信息分配因子α_i，并同时向局部滤波器S1和局部滤波器S2自适应重置过程噪声方差矩阵Q_i和估计误差协方差矩阵P_k-1,i；

所述信息分配因子

其中，(P_i,k-1)_nm表示k-1时刻估计误差协方差矩阵P_i,k-1的第n行、第m列的元素；

利用自适应重置后的Q_i和P_k-1,i实现微纳卫星多源信息姿态估计。

本发明实施例中自适应重置过程噪声方差矩阵Q_i和滤波过程估计误差协方差矩阵P_k-1,i为：

P_k-1,i＝α_i(P_k/k-1+I)

Q_i＝α_iQ

当k-1时刻局部滤波器估计误差协方差矩阵P_i,k-1初值选为零矩阵时，由于r_i的设计，使得r_i＞0，避免了信息分配因子α_i计算时的奇异问题。因此，本发明的改进型分布式自适应联邦卡尔曼滤波算法针对信息分配因子的修改，可有效避免当误差协方差矩阵初值选为零矩阵时，存在的奇异问题，同时保证了滤波的快速收敛及定姿精度。

本发明实施例一种可在轨全程应用的微纳卫星多源信息姿态确定方法，选择微型的姿态敏感器，如MEMS陀螺、磁强计和太阳敏感器；设计了一种改进型分布式自适应联邦卡尔曼滤波算法实现基于多源信息的微纳卫星的姿态确定，具体实施如下：

一、MEMS陀螺安装轴与微纳卫星惯性轴一致时，获取MEMS陀螺的测量值，也即微纳卫星本体坐标系b相对于惯性坐标系i的角速度ω_m，用于最终角速度的估计；

二、磁强计得到地磁场矢量的测量值，结合IGRF-12模型得到地磁矢量参考值，建立磁强计测量模型：具体实施如下

201利用国际地磁参考模型IGRF-12计算微纳卫星当前轨道位置下北东地地理坐标系下的地磁场矢量参考值，并经过地理坐标系、惯性坐标系及轨道坐标系之间的转换矩阵转换到轨道坐标系下的参考值B_o；

202磁强计安装轴与微纳卫星本体坐标系一致时，得到卫星本体坐标系下地磁场矢量测量值B_m，建立磁强计测量模型。

三、太阳敏感器得到测量坐标系下的太阳矢量测量值，结合太阳星历推导得到太阳矢量参考值，建立太阳敏感器的测量模型：具体实施如下

301微纳卫星上用的太阳敏感器主要是两轴敏感器。将两轴数字太阳敏感器的瞄准轴与微纳卫星本体坐标系-Z轴一致，得到两轴数字太阳敏感器的安装矩阵A_bm和太阳矢量在测量坐标系下的测量值S_m，经坐标转换获得本体坐标系下的太阳矢量测量值S_b；

302利用太阳星历计算惯性坐标系下的太阳矢量参考值，并经过坐标转换矩阵转换到轨道坐标系下，得到S_o，建立太阳敏感器的测量模型。

四、选择误差四元数Δq_bo的矢量部分Δq和陀螺角速率漂移估计误差Δb为状态量利用姿态运动学方程推导姿态确定系统的状态方程；

五、利用步骤二和步骤三中的磁强计和太阳敏感器的测量模型给出分布式自适应联邦卡尔曼滤波算法中的两个局部滤波器的观测信息；

六、利用改进型分布式自适应联邦卡尔曼滤波算法计算获得微纳卫星的姿态信息，实现基于多源信息的微纳卫星姿态解算与信息融合。改进型分布式自适应联邦卡尔曼滤波姿态确定系统主要包括局部滤波器S1根据反馈的姿态和自适应分配的协方差和过程噪声方差完成估计误差协方差的预测及基于磁强计的量测更新、局部滤波器S2根据反馈的姿态和自适应分配的协方差和过程噪声方差完成估计误差协方差的预测及基于太阳敏感器的量测更新；主滤波器实现姿态四元数的一步预测、协方差的一步预测、根据自适应分配因子进行信息分配、还有根据两个局部滤波器的输出信息进行信息融合、经过姿态和陀螺漂移修正后将最终姿态信息输出几大部分。

A.磁强计和太阳敏感器的测量模型

磁强计安装轴与微纳卫星本体坐标系一致时，得到卫星本体坐标系下地磁场矢量测量值B_m，利用国际地磁参考模型IGRF-12计算微纳卫星当前轨道位置下北东地地理坐标系下的地磁场矢量参考值，并经过转换矩阵转换到轨道坐标系下的参考值B_o；

建立磁强计测量模型为

B_m＝T_bo(q_bo)B_o+V_b

其中，V_b为磁强计的测量噪声；T_bo(q_bo)为卫星轨道系到本体系的方向余弦阵，由卫星当前姿态q_bo决定；B_o为轨道坐标系下的地磁场矢量。

本发明中选择两轴数字太阳敏感器作为姿态敏感器。将两轴数字太阳敏感器的瞄准轴与微纳卫星本体坐标系-Z轴一致，得到两轴数字太阳敏感器的安装矩阵A_bm和太阳矢量在测量坐标系下的测量值S_m，建立太阳敏感器的测量模型为

其中，V_s为太阳敏感器的测量噪声；T_bo(q_bo)为卫星轨道系到本体系的方向余弦阵，由卫星当前姿态决定；S_o为轨道坐标系下的太阳矢量。

给出S_o与本体坐标系下的太阳矢量S_b之间的关系为

S_b＝T_bo(q_bo)S_o

B.姿态确定系统的状态方程

给出四元数表示的姿态运动学方程如下

其中，q_bo＝[q₀ q₁ q₂ q₃]^T表示从轨道坐标系到卫星本体坐标系的四元数矢量。ω_bo＝[ω_x ω_y ω_z]^T。

卫星本体系相对于轨道坐标系的角速度矢量ω_bo可以表示为：

ω_bo＝ω_bI-T_bo(q_bo)ω_oI

其中，ω_bI是卫星本体坐标系相对于地心惯性坐标系的矢量，可以由MEMS陀螺测量得到；ω_oI＝[0 -ω_o 0]^T是轨道系相对于地心惯性系的角速度矢量；ω_o表示轨道运动角速度；表示为四元数乘法。

定义真实姿态q_bo与估计姿态间的误差四元数为：Δq_bo＝[Δq₀ Δq₁ Δq₂ Δq₃]^T＝[Δq₀ Δq^T]^T，由四元数合成法可得：

选择误差四元数Δq_bo的矢量部分Δq和陀螺角速率漂移估计误差Δb为状态量利用姿态运动学方程推导出姿态确定系统的状态方程为

其中：

将状态方程离散化后有

X_k＝Φ_k,k-1X_k-1+Γ_k,k-1W_k-1

其中：T为计算步长，k表示第k时刻，k-1表示第k-1时刻。

C.姿态确定系统的观测信息

利用磁强计测量模型与状态变量推导得到局部滤波器S1的观测信息为

Z_1k＝(B_b(k)-B_b(k/k-1))/2

其中：B_b(k)为k时刻卫星本体坐标系下地磁场矢量测量值；B_b(k/k-1)＝T_bo(q_bo(k/k-1))B_o(k)为根据IGRF数据和姿态估计值推算的本体坐标系下地磁场矢量；B_o(k)为k时刻卫星轨道坐标系下地磁场矢量；q_bo(k/k-1)为k-1时刻到k时刻的姿态四元数的预测值。

利用太阳敏感器测量模型与状态变量推导得到局部滤波器S2的观测信息为

Z_2k＝(S_b(k)-S_b(k/k-1))/2

其中：S_b(k)为k时刻卫星本体坐标系下太阳矢量测量值；S_b(k/k-1)＝T_bo(q_bo(k/k-1))S_o(k)为根据太阳星历和姿态估计值推算的本体坐标系下太阳矢量；S_o(k)为k时刻卫星轨道坐标系下太阳矢量。

如上给出的观测信息并非前、后两个采样时刻观测量的差值，而是当前时刻的观测量与根据当前姿态的预测值推算得到的观测量之间的差值。在观测信息的设置时，后者相比前者更准确，经过量测更新后所能获取的定姿精度更高。

D.改进型分布式联邦卡尔曼滤波过程

如图3所示，本发明中的改进型分布式自适应联邦卡尔曼滤波姿态确定系统主要包括局部滤波器S1完成估计误差协方差的预测及基于磁强计的量测更新、局部滤波器S2完成估计误差协方差的预测及基于太阳敏感器的量测更新；主滤波器实现姿态四元数的一步预测、协方差的一步预测、信息分配、信息融合、姿态和陀螺漂移修正及结果输出几个部分。具体实施步骤如下：

(1)根据估计的角速度和姿态运动学方程，利用RungeKutta方程估计姿态四元数

由k-1时刻的姿态估计四元数和估计角速度预测k时刻本体坐标系相对轨道坐标系的角速度如下

根据步骤四中的卫星姿态运动学方程，采用RungeKutta方程估计k时刻的姿态四元数为

其中：表示四元数乘法，为k-1时刻本体坐标系相对轨道坐标系的角速度，

(2)根据实施B中的离散状态信息，计算k时刻一步预测估计误差协方差矩阵P_k/k-1；

计算k时刻一步预测估计误差协方差矩阵P_k/k-1如下：

其中：P_k-1为k-1时刻的估计误差协方差矩阵；Q为主滤波器的过程噪声方差矩阵。

(3)计算信息分配因子α_i；

计算信息分配因子α_i如下

其中：

这种信息分配因子的计算，可有效避免当误差协方差矩阵初值选为零矩阵时，存在的奇异问题。

向局部滤波器S1和局部滤波器S2自适应重置过程噪声方差矩阵Q_i和估计误差协方差矩阵P_k-1,i。

P_k-1,i＝α_i(P_k/k-1+I)

Q_i＝α_iQ

(4)各局部滤波器计算k时刻的一步预测估计误差协方差矩阵P_i,k/k-1，并分别依据磁强计和太阳敏感器的观测信息进行量测更新；

其中：H_ik表示k时刻局部滤波器的观测矩阵；R_i表示局部滤波器的观测噪声方差矩阵。

(5)对局部滤波器S1和局部滤波器S2的状态估计和估计误差协方差进行数据融合，获得最终的状态估计和估计误差协方差矩阵P_k。

最终的状态估计和估计误差协方差矩阵P_k计算公式如下：

(6)利用状态估计修正估计姿态四元数和随机漂移估计值并结合MEMS陀螺的测量值ω_m得到最终的估计姿态四元数和估计角速度

实施例：

仿真实验：以某微纳卫星为例，该微纳卫星运行在轨道高度为546.8km的太阳同步轨道上，卫星转动惯量为diag[0.08845 0.1422 0.07518]kg·m²，降交点为10:30am，轨道倾角为97.62°。卫星在对地三轴稳定时，期望姿态为[0,0,0]°，MEMS陀螺的角度随机游走σ_υ＝0.03°/s^1/2，角速度随机游走σ_ζ＝0.0001°/s^3/2，磁强计测量噪声100nT，太阳敏感器的测量噪声为0.1°。利用本发明的微纳卫星姿态确定方法进行仿真实验。当微纳卫星对地三轴稳定转速率阻尼模式时，基于本发明改进型分布式自适应联邦卡尔曼滤波算法所确定的卫星的姿态信息如图4所示；基于本发明改进型分布式自适应联邦卡尔曼滤波算法所确定的卫星角速度信息如图5所示。当微纳卫星对地三轴稳定转太阳捕获巡航时，基于本发明改进型分布式自适应联邦卡尔曼滤波算法所确定的卫星姿态信息如图6所示；基于本发明改进型分布式自适应联邦卡尔曼滤波算法所确定的卫星角速度信息如图7所示。

由图4-图7可知，不管卫星处于哪种工作模式下或在某两种工作模式下切换，本发明的改进型分布式自适应联邦卡尔曼滤波算法所确定的卫星姿态和角速度都能够非常好的跟踪真实的姿态和角速度。仿真结果表明本发明具有很强的可实现性及较高的姿态估计精度。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种可在轨全程应用的微纳卫星多源信息姿态确定方法，其特征在于，包括如下步骤：

计算信息分配因子α_i，并同时向局部滤波器S1和局部滤波器S2自适应重置过程噪声方差矩阵Q_i和滤波过程估计误差协方差矩阵P_k-1,i，下角标i＝1,2分别对应两个局部滤波器；

所述信息分配因子

其中，(P_i,k-1)_nm表示k-1时刻估计误差协方差矩阵P_i,k-1的第n行、第m列的元素；

利用自适应重置后的Q_i和P_k-1,i实现微纳卫星多源信息姿态估计。

2.根据权利要求1所述可在轨全程应用的微纳卫星多源信息姿态确定方法，其特征在于，所述自适应重置过程噪声方差矩阵Q_i和滤波过程估计误差协方差矩阵P_k-1,i为：

P_k-1,i＝α_i(P_k/k-1+I)

Q_i＝α_iQ

其中，P_k/k-1表示k时刻主滤波器中一步预测估计误差协方差矩阵，Q为主滤波器的过程噪声方差矩阵。

3.根据权利要求1所述可在轨全程应用的微纳卫星多源信息姿态确定方法，其特征在于，所述局部滤波器S1的观测信息为：

Z_1k＝(B_b(k)-B_b(k/k-1))/2

B_b(k/k-1)＝T_bo(q_bo(k/k-1))B_o(k)

其中，B_b(k)为k时刻卫星本体坐标系下地磁场矢量测量值；B_o(k)为k时刻卫星轨道坐标系下地磁场矢量；q_bo(k/k-1)为k-1时刻到k时刻的姿态四元数的预测值；T_bo(·)为卫星轨道系到本体系的方向余弦阵。

4.根据权利要求1所述可在轨全程应用的微纳卫星多源信息姿态确定方法，其特征在于，所述局部滤波器S2的观测信息为：

Z_sk＝(S_b(k)-S_b(k/k-1))/2

S_b(k/k-1)＝T_bo(q_bo(k/k-1))S_o(k)

其中，S_b(k)为k时刻卫星本体坐标系下太阳矢量测量值；S_o(k)为k时刻卫星轨道坐标系下太阳矢量。