一种基于星光/地磁组合信息的小卫星自主导航系统及其导航方法

摘要

本发明提出一种基于星光/地磁组合信息的小卫星自主导航系统及其导航方法，属于组合导航技术领域。该自主导航系统包括星敏感器测量系统、地磁场测量系统、轨道确定系统和姿态确定系统；该导航方法包括：导航系统状态方程的建立、导航系统的量测方程建立和基于UKF的组合导航系统信息融合。本发明利用大视场星敏感器同时观测多颗导航星的星光矢量信息，弥补因磁强计测量噪声过大引起的缺陷，从而获得较高的导航精度。且本发明充分利用磁强计和星敏感器的测量值中包含的轨道和姿态信息，综合考虑各种因素的影响建立系统的量测方程，提高了滤波器对测量噪声的适应能力和稳定性，适用于对于采用小型廉价自主导航系统的低轨小卫星的导航。

说明书

技术领域

本发明涉及一种基于星光/地磁组合信息的小卫星自主导航方案设计系统及其导航方法， 属于组合导航技术领域。

背景技术

卫星自主导航是指卫星在不依赖地面站的条件下，仅依靠卫星上的测量设备实时地确定 卫星的位置、速度和姿态，在轨完成飞行任务所要求的功能或操作。当前卫星主要依靠地面 站进行跟踪与定轨，但由于小卫星组网中卫星数目较多，利用地面设备对卫星组网进行定位 存在着测量弧段较短和测量目标较多的困难，而且随着中低轨道小卫星的不断增加，传统的 地面测控设备难以满足实际应用需要。小卫星实现自主导航一方面有助于提高卫星自身的生 存能力；另一方面，还可以大大降低地面测控系统的成本。鉴于依赖地面站卫星导航的局限 性和卫星自主导航技术明显的优势，自主导航已成为卫星设计中的一个重要课题。

考虑到小卫星的特点，需要采用体积小、重量轻、低成本、可靠性高的导航系统实现其 轨道和姿态的确定。三轴磁强计具有自主性强、体积小、重量轻、成本低、性能可靠、没有 地域限制等优点，并能够提供全天侯、实时连续的导航信息。自Pasiaki等人提出利用磁强 计进行卫星自主导航以来，当前许多航天器都选用磁强计用于姿态的确定和控制。但是，由 于地磁场模型及磁强计测量的误差较大，目前单纯使用磁强计的地磁导航系统精度较低，因 此，需要利用其他的自主导航手段来弥补地磁导航的不足。

天文导航是利用天体敏感器对飞行器进行自主导航的一种方法，它具有隐蔽性好、定向 精度高、不受人工和自然电磁场的干扰、不向外辐射电磁波等优点，目前已经成功地应用在 各种航天器的自主导航系统中。星敏感器作为当前应用最广泛的天体敏感器，它是天文导航 系统的重要组成部分，它具有体积小、精度高、功耗低等特点，是理想的小型自主导航设备 之一，新型大视场星敏感器的出现，使得同时获取三颗及三颗以上导航星的星光矢量信息成 为可能，在不需要任何外部基准信息的前提下，即可实现高精度的自主导航，并且测量误差 不随时间积累。

目前地磁导航系统(GNS)能够提供的导航数据比较全面，不仅能够提供载体的位置信息， 而且能够提供载体的姿态信息，且数据输出率高，但地磁导航的精度较低，不能满足航天器 高精度的导航性能要求；天文导航(CNS)和地磁导航同属环境敏感导航系统，天文导航虽数 据率低，但其误差不随时间增加。因此，采用基于星光/地磁组合信息的导航方案(CNS/GNS) 可以满足小卫星自主导航的性能要求。

发明内容

针对现有技术中存在的问题，本发明提出一种基于星光/地磁组合信息的小卫星自主导航 系统及其导航方法，解决了磁强计测量精度较低而导致的地磁导航误差过大的问题，分析地 磁导航以及基于星敏感器的天文导航的特点和局限性，充分利用大视场星敏感器提供的多颗 导航星的星光矢量信息和三轴磁强计获得地球周围的地磁场信息，实现了二者的优势互补， 并通过UKF滤波方法进行信息融合，在提高精度的同时也提高了可靠性。

本发明提出一种基于星光/地磁组合信息的小卫星自主导航系统，包括星敏感器测量系 统、地磁场测量系统、轨道确定系统和姿态确定系统；

星敏感器测量系统包括大视场星敏感器、数据采集单元、星图识别单元、导航星历和 星光矢量获取单元；大视场星敏感器将测量得到的星光矢量信息提供给数据采集单元和星光 矢量获取单元；数据采集单元对来自大视场星敏感器的星光矢量信息进行数据预处理，将预 处理信息发送至星图识别单元；星图识别单元根据来自数据采集单元的预处理信息和导航星 历提供的参考信息进行星图匹配，获得星光矢量在惯性坐标系下的分量并将其提供给星光矢 量获取单元；星光矢量获取单元将大视场星敏感器测得星光矢量和星光矢量在惯性坐标系下 的分量提供给轨道确定系统和姿态确定系统；

地磁场测量系统包括三轴磁强计、国际地磁场模型和地磁矢量误差获取单元，三轴磁 强计通过测量地球表面以及近地空间的地磁场，获得当地地磁场矢量在卫星本体坐标系中的 分量，地磁场矢量的测量值；国际地磁场模型提供地磁场矢量的估计值；地磁矢量误差获取 单元利用三轴磁强计测得的地磁场矢量的测量值和国际地磁场模型得到的估计值计算二者的 地磁场矢量误差，并将其提供给轨道确定系统和姿态确定系统；

轨道确定系统包括卫星轨道动力学模型单元、β角获取单元、UKF滤波器A单元以及 卫星轨道信息获取单元；轨道确定系统根据轨道动力学模型单元建立用于轨道确定的状态方 程，并将该状态方程发送至UKF滤波器A单元；β角获取单元利用得到的地磁场矢量的测 量值和星光矢量信息计算二者之间的夹角，并将其作为观测量提供给UKF滤波器A单元； UKF滤波器A单元根据得到的观测量和用于轨道确定的状态方程，进行信息融合，实时估计 并校正卫星的轨道参数误差，并将获得的轨道参数提供给卫星轨道信息获取单元；卫星轨道 信息获取单元将得到的轨道参数提供给卫星，完成卫星的轨道确定；

姿态确定系统包括卫星姿态动力学模型单元、UKF滤波器B单元、姿态预估单元、姿态 误差估计单元和卫星姿态信息获取单元，姿态确定系统根据姿态动力学模型单元建立用于姿 态确定的状态方程；UKF滤波器B单元根据得到的测量信息和用于姿态确定的状态方程，进 行信息融合；姿态预估单元根据姿态动力学模型单元获得的姿态信息对卫星姿态进行预估计， 得到姿态预估信息；姿态误差估计单元根据来自UKF滤波器B单元的信息估计卫星的姿态 误差；姿态预估单元和姿态误差估计单元通过信息交流实现姿态校正，姿态预估单元将校正 后的姿态信息提供给卫星姿态信息获取单元；卫星姿态信息获取单元将得到的姿态信息提供 给卫星，实现卫星的姿态确定。

本发明提出一种基于星光/地磁组合信息的自主导航方法，具体包括以下步骤：

步骤一：导航系统状态方程的建立

定义导航系统的状态变量为

$\bar{X}={[\mathrm{\delta x},\mathrm{\delta y},\mathrm{\delta z},\delta v_x,\delta v_y,\delta v_z,\Delta \bar{w},\Delta \bar{q}]}^T---\left(1\right)$

式中，δx、δy、δz、δv_x、δv_y、δv_z分别表示卫星的位置和速度误差在惯性坐标系中的分量， 是轨道确定系统的状态量；为卫星本体坐标系相对于惯性空间的角速度误差；为误差 四元数矢量部分，和为姿态确定系统的状态量；

选取历元地心赤道惯性坐标系，考虑地球质心引力和J₂项非球形摄动的影响，将其它摄 动因素等效为高斯白噪声，卫星轨道动力学模型单元提供的轨道动力学模型为

$\left(\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x+{\omega }_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y+{\omega }_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z+{\omega }_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mu \frac{x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\omega }_{\mathrm{vx}}\\ \frac{{\mathrm{dv}}_y}{\mathrm{dt}}=-\mu \frac{y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\omega }_{\mathrm{vy}}\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mu \frac{z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\omega }_{\mathrm{vz}}\end{array}\right)---\left(2\right)$

其中，x、y、z分别表示卫星惯性坐标系中的位置；v_x、v_y、v_z分别表示卫星惯性坐标系 中的速度；为卫星到地心的距离；R_e为地球半径；μ为地心引力常数；J₂为 地球非球形二阶带谐引力摄动系数；ω_x、ω_y、ω_z分别为系统位置信息噪声；ω_vx、ω_vy、ω_vz分别为系统速度信息噪声；

定义卫星的轨道状态为公式(2)简化为其中， t表示时间；是轨道状态和时间t的函数；为系统噪声；

定义导航系统的位置和速度误差变量为

${\bar{x}}_o={\hat{\bar{X}}}_o-{\bar{X}}_o{[\mathrm{\delta x},\mathrm{\delta y},\mathrm{\delta z},\delta v_x,\delta v_y,\delta v_z]}^T---\left(3\right)$

其中，为的估计值；

则系统轨道确定部分的状态方程表示为

${\stackrel{·}{\stackrel{‾}{x}}}_o=f_o[{\bar{X}}_o,t]·{\bar{x}}_o+{\bar{\omega }}_o---\left(4\right)$

其中，为的微分；为相应的系统噪声；

采用四元数表示卫星姿态，则卫星姿态运动学模型单元提供的姿态动力学模型为

$\left(\begin{array}{c}\stackrel{·}{\stackrel{‾}{w}}={\bar{I}}^{-1}[{\bar{T}}_c+{\bar{T}}_d-\bar{w}\times (\bar{I}·\bar{w})\\ \stackrel{·}{q}=\frac{1}{2}q\otimes {\bar{w}}_{\mathrm{ob}}\\ {\bar{w}}_{\mathrm{ob}}=\bar{w}-C_o^b\left(q\right){\bar{w}}_{\mathrm{io}}\end{array}\right)---\left(5\right)$

其中，为卫星的控制力矩；为卫星的干扰力矩；为卫星的三轴转动惯量；为卫星 本体坐标系相对于轨道坐标系的角速度矢量；为卫星轨道角速度矢量，w_o为 卫星轨道角速率；为卫星相对于惯性空间的绝对角速度；q为卫星的姿态四元数；分 别为和q的微分；为四元数表示下的姿态矩阵；

根据四元数合成法则，并对公式(5)进行线性化，得到姿态确定系统的状态方程为

$\left(\begin{array}{c}\Delta \stackrel{·}{\stackrel{‾}{w}}={\bar{I}}^{-1}[\Delta {\bar{T}}_c+\Delta {\bar{T}}_d-\Delta \bar{w}\times (\bar{I}·\hat{\bar{w}})-\hat{\bar{w}}\times (\bar{I}·\Delta \bar{w})+2\Delta \bar{w}\times (\bar{I}·\Delta \bar{w})\\ \Delta \stackrel{·}{\stackrel{‾}{q}}=-[\hat{\bar{w}}\times ]·\Delta \bar{q}+\frac{1}{2}\Delta \bar{w}\end{array}\right)---\left(6\right)$

其中，表示卫星相对于惯性空间的绝对角速度，为其估计值；为卫星控制力矩的误 差；为卫星干扰力矩的误差；为卫星的三轴转动惯量；为卫星本体坐标系相对于惯 性空间的角速度误差，其微分为为误差四元数矢量部分，其微分为

根据公式(4)和(6)分别给出的轨道确定系统和姿态确定系统的状态方程，得到导航系统的 状态方程；

步骤二：导航系统的量测方程建立

分别选取大视场星敏感器测得的星光矢量信息、地磁矢量误差获取单元得到的地磁矢量 误差值以及三轴磁强计测得的地磁场矢量测量值与星光矢量信息的夹角作为系统观测量，建 立量测方程；

a、基于星光矢量信息的量测方程建立

星敏感器测量系统的星光矢量获取单元同时输出多颗导航星的星光矢量信息，姿态确定 系统利用来自星光矢量获取单元的两颗导航恒星的星光矢量信息L_m1、L_m2进行卫星的姿态确 定；大视场星敏感器测得的星光矢量信息L_m可表示为

$L_m=C_o^b\left(q\right)C_i^o{L}_i+{\upsilon }_L---\left(7\right)$

式中，L_i为星光矢量在惯性系中的分量，利用星图识别单元和导航星历获取；υ_L为星敏感器 测量误差；为四元数表示下的姿态矩阵；为惯性系到轨道坐标系的转换矩阵；q表 示卫星的姿态四元数；结合公式(7)表示为

$L_m={\left[C_o^b\left(\mathrm{\Delta q}\right)\right]}^T{C}_o^b\left(\hat{q}\right)C_i^o{L}_i+{\upsilon }_L$

(8)

$=(I_{3\times 3}+2[\Delta \bar{q}\times \left]\right)C_o^b\left(\hat{q}\right)C_i^o{L}_i+{\upsilon }_L={\hat{L}}_b+2[{\hat{L}}_b\times ]\Delta \bar{q}+{\upsilon }_L$

式中，为星光矢量在本体坐标系中的估计值；为估计姿态四元数下的姿态矩阵； 为误差四元数下的姿态矩阵；I_3×3表示单位矩阵；为误差四元数的矢量部分；

得到基于星光矢量的量测方程Z₁为

$Z_1=\left(\begin{array}{c}{\hat{L}}_{b1}+2[{\hat{L}}_{b1}\times ]\Delta \bar{q}+{\upsilon }_{L1}\\ {\hat{L}}_{b2}+2[{\hat{L}}_{b2}\times ]\Delta \bar{q}+{\upsilon }_{L2}\end{array}\right)---\left(9\right)$

其中，为导航星1的星光矢量在本体坐标系中的估计值；υ_L1为导航星1的星光矢量的测 量误差；为导航星2的星光矢量在本体坐标系中的估计值；υ_L2为导航星2的星光矢量的 测量误差；

b、基于地磁矢量误差的量测方程建立

三轴磁强计输出的地磁场在卫星本体系下的分量B_m可表示为

$B_m=C_i^b{C}_n^i{B}_n+{\upsilon }_{B1}---\left(10\right)$

其中，B_n为地磁场矢量在地理坐标系中的真实值；分别为地理系到惯性系的转换矩 阵、惯性系到本体系的转换矩阵；υ_B1为三轴磁强计的地磁场矢量测量误差；

地磁场矢量的估计值在本体系中的分量为

${\hat{B}}_b={\hat{C}}_i^b{\hat{C}}_e^i{\hat{B}}_n+{\upsilon }_{B2}---\left(11\right)$

其中，为B_n的估计值，根据国际参考地磁场(IGRF)模型求得；分别为和的 估计值；υ_B2为地磁场模型误差；

则地磁矢量误差获取单元获得的观测量ΔB为

$\mathrm{\Delta B}={\hat{B}}_b-B_m={\hat{C}}_i^b{\hat{C}}_n^i{\hat{B}}_n-C_i^b{C}_n^i{B}_n+{\upsilon }_{B2}-{\upsilon }_{B1}---\left(12\right)$

定义$\Delta \left(C_i^b{C}_n^i{B}_n\right)={\hat{C}}_i^b{\hat{C}}_n^i{\hat{B}}_n-C_i^b{C}_n^i{B}_n,$υ_B＝υ_B2-υ_B1，则公式(12)简化为

$\mathrm{\Delta B}=\Delta \left(C_i^b{C}_n^i{B}_n\right)+{\upsilon }_B---\left(13\right)$

根据微分运算法则，公式(13)右端第一项表示为

$\Delta \left(C_i^b{C}_n^i{B}_n\right)=C_i^b\Delta \left(C_n^i{B}_n\right)+\Delta C_i^b\left(C_n^i{B}_n\right)---\left(14\right)$

根据地磁场矢量在地理坐标系中的分量B_n的定义，建立同轨道误差状态的关 系表达式H_o的表达式为

$H_o[{\bar{X}}_o,t]=C_i^b·\frac{\partial \left[C_n^i(\theta ,\lambda )B_n(r,\theta ,\lambda )\right]}{\partial {\bar{X}}_o}---\left(15\right)$

其中，

$\frac{\partial \left[C_n^i(\theta ,\lambda )B_n(r,\theta ,\lambda )\right]}{\partial {\bar{X}}_o}$

$=(\frac{\partial C_n^i}{\partial \theta }{B}_n+C_n^i\frac{\partial B_n}{\partial \theta })\frac{\partial \theta }{\partial {\bar{X}}_o}+(\frac{\partial C_n^i}{\partial \lambda }{B}_n+C_n^i\frac{\partial B_n}{\partial \lambda })\frac{\partial \lambda }{\partial {\bar{X}}_o}+C_n^i(\frac{\partial B_n}{\partial {\bar{X}}_o}+C_n^i\frac{\partial B_n}{\partial r}\frac{\partial r}{\partial {\bar{X}}_o})$

公式(14)右端第二项改写为

$\Delta C_i^b\left(C_n^i{B}_n\right)=\Delta C_i^b{B}_i---\left(16\right)$

其中，B_i为真实的地磁场矢量在惯性系中的分量；是由卫星的姿态误差引起，其表达式 为

$\Delta C_i^b={\hat{C}}_i^b-C_i^b---\left(17\right)$

根据姿态四元数为小量，可以得到

$\Delta C_i^b={\hat{C}}_i^b-C_i^b=(I_{3\times 3}-2[\Delta \bar{q}\times \left]\right)C_i^b-C_i^b=-2[\Delta \bar{q}\times ]C_i^b---\left(18\right)$

其中，I_3×3表示单位矩阵；为误差四元数的矢量部分；

将式(18)带入式(16得到

$\Delta C_i^b{B}_i=(-2[\Delta \bar{q}\times \left]C_i^b\right)B_i=-2[\Delta \bar{q}\times ]B_b=2[B_b\times ]\Delta \bar{q}---\left(19\right)$

其中，B_b为真实的地磁场矢量在体坐标系中的分量；

根据公式(13)～(19)，公式(12)表示为

$\mathrm{\Delta B}=2[B_b\times ]\Delta \bar{q}+H_o{\bar{x}}_o+{\upsilon }_B---\left(20\right)$

由于B_b为未知量，采用三轴磁强计的测量值B_m代替；则导航系统的观测方程Z₂为

$Z_2=2[B_m\times ]\Delta \bar{q}+H_o{\bar{x}}_o+{\upsilon }_B---\left(21\right)$

其中，H_o表示轨道误差状态的量测矩阵；υ_B为相应的量测噪声；

c、基于星光矢量与地磁矢量夹角的量测方程建立

β角获取单元选取磁强计测得的地磁矢量B_m与星敏感器测量得到的星光矢量L_m3之间的 夹角作为观测量，相应的量测方程Z₃为

$Z_3=\beta =\arccos \left(\frac{L_{m3}^T{B}_m}{\left|\right|B_m\left|\right|}\right)+{\upsilon }_{\beta }---\left(22\right)$

其中，||B_m||为B_m的模；υ_β夹角的测量误差；

根据公式(9)、(21)以及(22)得到导航系统的量测方程；

步骤三：基于UKF的组合导航系统信息融合

UKF滤波器A单元利用轨道动力学模型单元建立轨道确定系统的状态方程，将来自地磁 矢量误差获取单元的地磁矢量误差ΔB和β角获取单元获得的夹角β作为观测量，在UKF滤 波的框架下对卫星的轨道误差状态进行实时估计，并将其发送给卫星轨道信息获取单元，对 卫星的轨道参数进行校正，实现轨道确定；

UKF滤波器B单元利用姿态动力学模型单元建立姿态确定系统的状态方程，将来自地磁 矢量误差获取单元的地磁矢量误差ΔB和星光矢量获取单元的星光矢量信息L_m1和L_m2作为观 测量，在UKF滤波的框架下对卫星的姿态误差状态进行估计，并将其提供给姿态误差估计单 元；姿态预估单元利用来自姿态动力学模型单元的信息进行姿态预估，并通过信息交换与姿 态误差估计单元一起进行姿态校正，把校正后姿态信息输入到姿态信息获取单元中，从而完 成卫星的姿态确定。

本发明的优点在于：

(1)本发明提出一种基于星光/地磁组合信息的小卫星自主导航系统及其导航方法，利用 大视场星敏感器同时观测多颗导航星的星光矢量信息，弥补因磁强计测量噪声过大引起的缺 陷，从而获得较高的导航精度；

(2)本发明提出一种基于星光/地磁组合信息的小卫星自主导航系统及其导航方法，充分 利用磁强计和星敏感器的测量值中包含的轨道和姿态信息，综合考虑各种因素的影响建立系 统的量测方程，提高了滤波器对测量噪声的适应能力和稳定性；

(3)本发明提出一种基于星光/地磁组合信息的小卫星自主导航系统及其导航方法，根据 卫星的轨道和姿态动力学模型，推导并建立了能够同时完成轨道和姿态确定的系统状态方程， 并通过UKF算法进行信息融合完成轨道和姿态的确定；

(4)本发明提出一种基于星光/地磁组合信息的小卫星自主导航系统及其导航方法，具有 结构简单、体积小、系统可靠性高等特点，能够满足低轨小卫星自主导航要求，对于采用小 型廉价自主导航系统的低轨小卫星而言，是一种非常有发展前景的卫星自主导航方法。

附图说明

图1为本发明提出的一种基于星光/地磁组合信息的小卫星自主导航系统结构示意图；

图2为本发明中地磁场球面坐标系及各地磁要素之间的关系图；

图中：

1-星敏感器测量系统    2-地磁场测量系统   3-轨道确定系统    4-姿态确定系统

101-大视场星敏感器    201-三轴磁强计     301-轨道动力学    401-姿态动力学模

                                         模型单元          型单元

102-数据采集单元      202-国际地磁场     302-β角获取单    402-UKF滤波器B

                      (IGRF)模型         元                单元

103-星图识别单元      203-地磁矢量误差   303-UKF滤波器     403-姿态预估单元

                      获取单元           A单元

                                         304-卫星轨道信    404-姿态误差估计

104-导航星历

                                         息获取单元        单元

105-星光矢量获取单元                                       405-卫星姿态信息

                                                         获取单元

具体实施方式

下面将结合附图对本发明作进一步的详细说明。

本发明提出一种基于星光/地磁组合信息的小卫星自主导航系统及其导航方法，如图1 所示，该自主导航系统的设计分为轨道确定和姿态确定两部分。根据卫星的轨道动力学模型 和姿态动力学模型建立系统状态方程；以星敏感器测得的星光矢量(L_m1、L_m2和L_m3)、星光矢 量与磁强计测得的地磁矢量(B_m)之间的夹角(β)、地磁矢量的测量值(B_m)与估计值之差 (ΔB)作为系统观测量；利用UKF滤波算法估计卫星的轨道和姿态参数，实现小卫星的自主 导航。

本发明提出的一种基于星光/地磁组合信息的小卫星自主导航系统，包括星敏感器测量系 统1、地磁场测量系统2、轨道确定系统3和姿态确定系统4。

星敏感器测量系统1主要包括大视场星敏感器101、数据采集单元102、星图识别单 元103、导航星历104和星光矢量获取单元105。大视场星敏感器101能在同一时刻观测 得到三颗以上导航恒星的星光矢量信息，并将测量得到的星光矢量信息提供给数据采集单元 102和星光矢量获取单元105；数据采集单元102对来自大视场星敏感器101的星光矢量 信息进行数据预处理，将预处理信息发送至星图识别单元103；星图识别单元103根据来自 数据采集单元102的预处理信息和导航星历104提供的参考信息进行星图匹配，获得星光 矢量在惯性坐标系下的分量并将其提供给星光矢量获取单元105；星光矢量获取单元105将 大视场星敏感器101测得星光矢量和星光矢量在惯性坐标系下的分量提供给轨道确定系统3 和姿态确定系统4。

地磁场测量系统2主要包括三轴磁强计201、国际地磁场(IGRF)模型202和地磁矢量 误差获取单元203。在轨道确定和姿态确定过程中，三轴磁强计201通过测量地球表面以及 近地空间的地磁场，获得当地地磁场矢量在卫星本体坐标系中的分量，即地磁场矢量的测量 值；国际地磁场(IGRF)模型202可以提供地磁场矢量的估计值；地磁矢量误差获取单元 203利用三轴磁强计201测得的地磁场矢量的测量值和国际地磁场(IGRF)模型202得到 的估计值计算二者的地磁场矢量误差，并将其提供给轨道确定系统3和姿态确定系统4。

轨道确定系统3包括卫星轨道动力学模型单元301、β角获取单元302、UKF滤波器 A单元303以及卫星轨道信息获取单元304。在轨道确定过程中，轨道确定系统3根据轨道 动力学模型单元301建立用于轨道确定的状态方程，并将该状态方程发送至UKF滤波器A 单元303；β角获取单元302利用得到的地磁场矢量的测量值和星光矢量信息计算二者之间 的夹角，并将其作为观测量提供给UKF滤波器A单元303；UKF滤波器A单元303根据 得到的观测量和用于轨道确定的状态方程，进行信息融合，实时估计并校正卫星的轨道参数 误差，并将获得的轨道参数提供给卫星轨道信息获取单元304；卫星轨道信息获取单元304 将得到的轨道参数提供给卫星，完成卫星的轨道确定。

姿态确定系统4主要包括卫星姿态动力学模型单元401、UKF滤波器B单元402、姿 态预估单元403、姿态误差估计单元404和卫星姿态信息获取单元405。在姿态确定过程中， 姿态确定系统4根据姿态动力学模型单元401建立用于姿态确定的状态方程；UKF滤波器 B单元402根据得到的测量信息和用于姿态确定的状态方程，进行信息融合；姿态预估单元 403根据姿态动力学模型单元401获得的姿态信息对卫星姿态进行预估计，得到姿态预估信 息；姿态误差估计单元404根据来自UKF滤波器B单元402的信息估计卫星的姿态误差； 姿态预估单元403和姿态误差估计单元404通过信息交流实现姿态校正，姿态预估单元403 将校正后的姿态信息提供给卫星姿态信息获取单元405；卫星姿态信息获取单元405将得到 的姿态信息提供给卫星，实现卫星的姿态确定。

本发明提出一种基于星光/地磁组合信息的自主导航方法，具体包括以下步骤：

步骤一：导航系统状态方程的建立

定义导航系统的状态变量为

$\bar{X}={[\mathrm{\delta x},\mathrm{\delta y},\mathrm{\delta z},\delta v_x,\delta v_y,\delta v_z,\Delta \bar{w},\Delta \bar{q}]}^T---\left(1\right)$

式中，δx、δy、δz、δv_x、δv_y、δv_z分别表示卫星的位置和速度误差在惯性坐标系中的分量， 是轨道确定系统3的状态量；为卫星本体坐标系相对于惯性空间的角速度误差；为误 差四元数矢量部分，和为姿态确定系统4的状态量。

选取历元(J2000.0)地心赤道惯性坐标系，主要考虑地球质心引力和J₂项非球形摄动的 影响，将其它摄动因素(指除了地球质心引力和J₂项非球形摄动以外的摄动)等效为高斯白 噪声，则卫星轨道动力学模型单元301提供的轨道动力学模型为

$\left(\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x+{\omega }_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y+{\omega }_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z+{\omega }_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mu \frac{x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\omega }_{\mathrm{vx}}\\ \frac{{\mathrm{dv}}_y}{\mathrm{dt}}=-\mu \frac{y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\omega }_{\mathrm{vy}}\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mu \frac{z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\omega }_{\mathrm{vz}}\end{array}\right)---\left(2\right)$

式中，x、y、z分别表示卫星惯性坐标系中的位置；v_x、v_y、v_z分别表示卫星惯性坐标系 中的速度；为卫星到地心的距离；R_e为地球半径；μ为地心引力常数；J₂为 地球非球形二阶带谐引力摄动系数；ω_x、ω_y、ω_z分别为系统位置信息噪声；ω_vx、ω_vy、ω_vz分别为系统速度信息噪声。

定义卫星的轨道状态为公式(2)可以简写为其中，t表示时间；是轨道状态和时间t的函数；为系统噪声。

定义导航系统的位置和速度误差变量为

${\bar{x}}_o={\hat{\bar{X}}}_o-{\bar{X}}_o{[\mathrm{\delta x},\mathrm{\delta y},\mathrm{\delta z},\delta v_x,\delta v_y,\delta v_z]}^T---\left(3\right)$

式中，为的估计值。

则系统轨道确定部分3的状态方程可表示为

${\stackrel{·}{\stackrel{‾}{x}}}_o=f_o[{\bar{X}}_o,t]·{\bar{x}}_o+{\bar{\omega }}_o---\left(4\right)$

式中，为的微分；为相应的系统噪声。

采用四元数表示卫星姿态，则卫星姿态运动学模型单元401提供的姿态动力学模型为

$\left(\begin{array}{c}\stackrel{·}{\stackrel{‾}{w}}={\bar{I}}^{-1}[{\bar{T}}_c+{\bar{T}}_d-\bar{w}\times (\bar{I}·\bar{w})\\ \stackrel{·}{q}=\frac{1}{2}q\otimes {\bar{w}}_{\mathrm{ob}}\\ {\bar{w}}_{\mathrm{ob}}=\bar{w}-C_o^b\left(q\right){\bar{w}}_{\mathrm{io}}\end{array}\right)---\left(5\right)$

式中，为卫星的控制力矩；为卫星的干扰力矩；为卫星的三轴转动惯量；为卫星 本体坐标系相对于轨道坐标系的角速度矢量；为卫星轨道角速度矢量，w_o为 卫星轨道角速率；为卫星相对于惯性空间的绝对角速度；q为卫星的姿态四元数；分 别为和q的微分；为四元数表示下的姿态矩阵。

根据四元数合成法则，并对表达式(5)进行线性化可以得到姿态确定系统4的状态方程为

$\left(\begin{array}{c}\Delta \stackrel{·}{\stackrel{‾}{w}}={\bar{I}}^{-1}[\Delta {\bar{T}}_c+\Delta {\bar{T}}_d-\Delta \bar{w}\times (\bar{I}·\hat{\bar{w}})-\hat{\bar{w}}\times (\bar{I}·\Delta \bar{w})+2\Delta \bar{w}\times (\bar{I}·\Delta \bar{w})\\ \Delta \stackrel{·}{\stackrel{‾}{q}}=-[\hat{\bar{w}}\times ]·\Delta \bar{q}+\frac{1}{2}\Delta \bar{w}\end{array}\right)---\left(6\right)$

其中，表示卫星相对于惯性空间的绝对角速度，为其估计值；为卫星控制力矩的误 差；为卫星干扰力矩的误差；为卫星的三轴转动惯量；为卫星本体坐标系相对于惯 性空间的角速度误差，其微分为为误差四元数矢量部分，其微分为

由上式可知，误差四元数自然降阶为四元数矢量部分，从而避免了由于四元数正交性约 束所造成的奇异性问题。卫星的自主导航分为设计分为轨道确定和姿态确定两部分，根据式 (4)和(6)分别给出的轨道确定系统和姿态确定系统的状态方程，即可得到导航系统的状态方 程。

步骤二：导航系统的量测方程建立

分别选取大视场星敏感器101测得的星光矢量信息、地磁矢量误差获取单元203得到 的地磁矢量误差值以及三轴磁强计201测得的地磁场矢量测量值与星光矢量信息的夹角作为 系统观测量，建立量测方程。

a、基于星光矢量信息的量测方程建立

根据图1可知，星敏感器测量系统1的星光矢量获取单元105可以同时输出多颗导航星 的星光矢量信息。姿态确定系统4利用来自星光矢量获取单元105的两颗导航恒星(导航星1 和导航星2)的星光矢量信息(L_m1、L_m2)进行卫星的姿态确定。大视场星敏感器101测得的星 光矢量信息L_m可表示为

$L_m=C_o^b\left(q\right)C_i^o{L}_i+{\upsilon }_L---\left(7\right)$

式中，L_i为星光矢量在惯性系中的分量，可利用星图识别单元103和导航星历104获取；υ_L为星敏感器测量误差；为四元数表示下的姿态矩阵；为惯性系到轨道坐标系的转换 矩阵；q表示卫星的姿态四元数；结合式(7)可表示为

$L_m={\left[C_o^b\left(\mathrm{\Delta q}\right)\right]}^T{C}_o^b\left(\hat{q}\right)C_i^o{L}_i+{\upsilon }_L$

(8)

$=(I_{3\times 3}+2[\Delta \bar{q}\times \left]\right)C_o^b\left(\hat{q}\right)C_i^o{L}_i+{\upsilon }_L={\hat{L}}_b+2[{\hat{L}}_b\times ]\Delta \bar{q}+{\upsilon }_L$

式中，为星光矢量在本体坐标系中的估计值；为估计姿态四元数下的姿态矩阵； 为误差四元数下的姿态矩阵；I_3×3表示单位矩阵；为误差四元数的矢量部分。

得到基于星光矢量的量测方程Z₁为

$Z_1=\left(\begin{array}{c}{\hat{L}}_{b1}+2[{\hat{L}}_{b1}\times ]\Delta \bar{q}+{\upsilon }_{L1}\\ {\hat{L}}_{b2}+2[{\hat{L}}_{b2}\times ]\Delta \bar{q}+{\upsilon }_{L2}\end{array}\right)---\left(9\right)$

其中，为导航星1的星光矢量在本体坐标系中的估计值；υ_L1为导航星1的星光矢量的测 量误差；为导航星2的星光矢量在本体坐标系中的估计值；υ_L2为导航星2的星光矢量的 测量误差。

b、基于地磁矢量误差的量测方程建立

三轴磁强计201输出的地磁场在卫星本体系下的分量B_m可表示为

$B_m=C_i^b{C}_n^i{B}_n+{\upsilon }_{B1}---\left(10\right)$

式中，B_n为地磁场矢量在地理坐标系中的真实值；分别为地理系到惯性系的转换矩 阵、惯性系到本体系的转换矩阵；υ_B1为三轴磁强计201的地磁场矢量测量误差。

地磁场矢量的估计值在本体系中的分量为

${\hat{B}}_b={\hat{C}}_i^b{\hat{C}}_e^i{\hat{B}}_n+{\upsilon }_{B2}---\left(11\right)$

其中，为B_n的估计值，可根据国际参考地磁场(IGRF)模型202求得；分别为和 的估计值；υ_B2为地磁场模型误差。

则地磁矢量误差获取单元203获得的观测量ΔB为

$\mathrm{\Delta B}={\hat{B}}_b-B_m={\hat{C}}_i^b{\hat{C}}_n^i{\hat{B}}_n-C_i^b{C}_n^i{B}_n+{\upsilon }_{B2}-{\upsilon }_{B1}---\left(12\right)$

定义$\Delta \left(C_i^b{C}_n^i{B}_n\right)={\hat{C}}_i^b{\hat{C}}_n^i{\hat{B}}_n-C_i^b{C}_n^i{B}_n,$υ_B＝υ_B2-υ_B1，则式(12)可简化为

$\mathrm{\Delta B}=\Delta \left(C_i^b{C}_n^i{B}_n\right)+{\upsilon }_B---\left(13\right)$

根据微分运算法则，式(13)右端第一项可以表示为

$\Delta \left(C_i^b{C}_n^i{B}_n\right)=C_i^b\Delta \left(C_n^i{B}_n\right)+\Delta C_i^b\left(C_n^i{B}_n\right)---\left(14\right)$

图2给出了地磁场球面坐标系及各地磁要素之间的关系图。图中，O_iX_iY_iZ_i表示地心惯 性系；O_eX_eY_eZ_e表示地球坐标系；G为Greenwich子午线赤经；B_r、B_θ、B_λ为地磁场矢 量在地磁场球面坐标系三个轴上的分量；r为卫星距地心的距离，λ为地理经度，θ为地心余 纬；ONET表示当地地理坐标系，O为卫星的质心，ON指向地理北，OE指向地理东，OT垂 直于ONE平面且指向地心。B_n为地磁场矢量在地理坐标系中的分量，B_H为B_n在ONE平面 内的分量，称为水平分量；X为B_H在ON轴上的分量，称为北向分量；Y为B_H在OE轴上的 分量，称为东向分量；Z为B_n在OT轴上的分量，称为垂直强度或垂直分量；D为B_H偏离地 理北的角度，称为磁偏角，向东偏为正；I为B_n与水平面ONE之间的夹角，称为磁倾角，向 下倾为正。

根据图2给出的地磁场矢量在地理坐标系中的分量B_n的定义可知，地理系到惯性系的转 换矩阵B_n与卫星的位置信息(卫星的地心距r、地理经度λ和地心余纬θ)有关，即等式 (14)右端第一项主要由轨道参数误差引起。因此，可以建立同轨道误差状态的 关系表达式H_o的表达式为

$H_o[{\bar{X}}_o,t]=C_i^b·\frac{\partial \left[C_n^i(\theta ,\lambda )B_n(r,\theta ,\lambda )\right]}{\partial {\bar{X}}_o}---\left(15\right)$

其中，

$\frac{\partial \left[C_n^i(\theta ,\lambda )B_n(r,\theta ,\lambda )\right]}{\partial {\bar{X}}_o}$

$=(\frac{\partial C_n^i}{\partial \theta }{B}_n+C_n^i\frac{\partial B_n}{\partial \theta })\frac{\partial \theta }{\partial {\bar{X}}_o}+(\frac{\partial C_n^i}{\partial \lambda }{B}_n+C_n^i\frac{\partial B_n}{\partial \lambda })\frac{\partial \lambda }{\partial {\bar{X}}_o}+C_n^i(\frac{\partial B_n}{\partial {\bar{X}}_o}+C_n^i\frac{\partial B_n}{\partial r}\frac{\partial r}{\partial {\bar{X}}_o})$

式(14)右端第二项与卫星的姿态测量矩阵有关，其可改写为

$\Delta C_i^b\left(C_n^i{B}_n\right)=\Delta C_i^b{B}_i---\left(16\right)$

其中，B_i为真实的地磁场矢量在惯性系中的分量；主要是由卫星的姿态误差引起，其表 达式为

$\Delta C_i^b={\hat{C}}_i^b-C_i^b---\left(17\right)$

根据姿态四元数为小量，可以得到

$\Delta C_i^b={\hat{C}}_i^b-C_i^b=(I_{3\times 3}-2[\Delta \bar{q}\times \left]\right)C_i^b-C_i^b=-2[\Delta \bar{q}\times ]C_i^b---\left(18\right)$

其中，I_3×3表示单位矩阵；为误差四元数的矢量部分。

将式(18)带入式(16)，可得

$\Delta C_i^b{B}_i=(-2[\Delta \bar{q}\times \left]C_i^b\right)B_i=-2[\Delta \bar{q}\times ]B_b=2[B_b\times ]\Delta \bar{q}---\left(19\right)$

其中，B_b为真实的地磁场矢量在体坐标系中的分量。

根据式(13)～(19)，式(12)可表示为

$\mathrm{\Delta B}=2[B_b\times ]\Delta \bar{q}+H_o{\bar{x}}_o+{\upsilon }_B---\left(20\right)$

由于B_b为未知量，在实际过程中可采用三轴磁强计201的测量值B_m代替。则导航系统 的观测方程Z₂为

$Z_2=2[B_m\times ]\Delta \bar{q}+H_o{\bar{x}}_o+{\upsilon }_B---\left(21\right)$

其中，H_o表示轨道误差状态的量测矩阵；υ_B为相应的量测噪声。

c、基于星光矢量与地磁矢量夹角的量测方程建立

β角获取单元302选取磁强计201测得的地磁矢量B_m与星敏感器101测量得到的星光 矢量L_m3之间的夹角作为观测量，相应的量测方程Z₃为

$Z_3=\beta =\arccos \left(\frac{L_{m3}^T{B}_m}{\left|\right|B_m\left|\right|}\right)+{\upsilon }_{\beta }---\left(22\right)$

其中，||B_m||为B_m的模；υ_β夹角的测量误差。

根据式(9)、(21)以及(22)，即可得到导航系统的量测方程。

步骤三：基于UKF的组合导航系统信息融合

在基于星光/地磁组合信息的自主导航系统中，三轴磁强计201和大视场星敏感器101 的测量信息均与卫星的轨道和姿态信息有关，因此导航系统可以利用UKF滤波算法同时实现 卫星的轨道确定和姿态确定。

UKF滤波器A单元303利用轨道动力学模型单元301建立轨道确定系统的状态方程， 将来自地磁矢量误差获取单元203的地磁矢量误差ΔB和β角获取单元302获得的夹角β作 为观测量，在UKF滤波的框架下对卫星的轨道误差状态进行实时估计，并将其发送给卫星轨 道信息获取单元304，对卫星的轨道参数进行校正，实现轨道确定。

UKF滤波器B单元402利用姿态动力学模型单元401建立姿态确定系统的状态方程， 将来自地磁矢量误差获取单元203的地磁矢量误差ΔB和星光矢量获取单元105的星光矢量 信息(L_m1、L_m2)作为观测量，在UKF滤波的框架下对卫星的姿态误差状态进行估计，并将其 提供给姿态误差估计单元404；姿态预估单元403利用来自姿态动力学模型单元401的信 息进行姿态预估，并通过信息交换与姿态误差估计单元404一起进行姿态校正，把校正后姿 态信息输入到姿态信息获取单元405中，从而完成卫星的姿态确定。