微小卫星磁传感器多误差模型在轨迭代标定方法

摘要

本发明公开了一种基于多参数的二步微小卫星磁传感器多误差模型在轨迭代标定方法。首先在光照区利用太阳敏感器和陀螺计算得到卫星俯仰角，然后以地面标定或上一次测量的标度因数为标准，只标定零偏误差，最后重新建立多误差标定模型，采用最佳拟合的方法重新标定出零偏和标度因数误差。本发明通过标定磁强计的零偏和标度因数误差，并在姿态解算中予以补偿，从而有效地提高卫星导航精度；其不需要专门针对标定算法设计新的滤波器，也不需要实时计算，有效地降低了卫星运行成本和负载。

说明书

技术领域

本发明涉及一种微小卫星磁传感器多误差模型在轨迭代标定方法，尤其是一 种基于多参数的二步微小卫星磁传感器多误差模型在轨迭代标定方法。

背景技术

磁强计是微小卫星姿态测量的主要部件，磁强计的测量精度是影响微小卫星 定姿性能的重要指标。

卫星的剩磁干扰是星上磁强计的主要误差来源，控制去除剩磁干扰的方法主 要包括：(1)、磁补法。即将发射前的剩磁测量与补偿的方法相结合对剩磁量加 以控制，但该方法作用有限，经过磁补以后的剩磁量一般还有900nT左右，对姿 态测量将产生较大的误差；(2)、标定法。就是确定磁强计的误差参数，以便在 姿态解算中用于计算补偿。传统的在轨标定算法中多数只标定零偏误差，但标度 因数误差也应考虑，还有学者提出用多阶滤波器进行全误差模型标定，需要实时 计算，增加了系统的计算负担。

2011年12月，中国惯性技术学报,第19卷第6期，614-617页，作者：华冰， 郁丰，公开了一种基于光照区太阳敏感器和陀螺辅助修正的微小卫星磁测技术， 其采用太阳敏感器和陀螺对磁强计误差进行辅助测量与修正，推导了磁强计误差 估计方法，在光照区以太阳敏感器与陀螺输出作为俯仰滤波器观测量，估计出卫 星俯仰角度和角速度，再采用最小二乘方法，利用滤波输出量对磁强计误差进行 估计，估计的结果进入滤波器对磁场输出进行测量修正。该方法，姿态角精度提 高了1°左右，角速度精度最高提高了0.003(°)/s左右，并增强了卫星稳定性。 但是该方法只能标定零偏误差，无法标定标度因数误差，导致只能标定简化的磁 强计模型，实用精度并不高。

发明内容

本发明所要解决的技术的问题在于克服现有技术的缺陷，提供一种通过同时 标定零偏误差和标度因数误差来提高导航定位精度的微小卫星磁传感器多误差 模型在轨迭代标定方法。

为了解决上述技术问题，本发明提供的微小卫星磁传感器多误差模型在轨 迭代标定方法，包括以下步骤：

1)、在光照区利用太阳敏感器和陀螺计算得到卫星俯仰角θ；

2)、以地面标定或者上一次测量的标度因数为标准，建立零偏误差粗标定 模型：

H_iE＝G_i      (14)

式(14)中,H_i为量测量和已知量计算得到的用于标定的矢量， H_i＝[k_z(b_xosinθ+b_zocosθ)k_x(b_zosinθ-b_xocosθ)]；E为磁强计零偏误差向 量，G_i为量测量和已知量计算得到的用于标定的标量， $G_i=k_z{\tilde{b}}_x(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )+k_x{\tilde{b}}_z(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta );$其中，θ为卫星俯 仰角，为X轴磁强计测量值，为X轴磁强计零偏误差，k_x为X轴磁强计标度 因数，为Z轴磁强计测量值，为Z轴磁强计零偏误差，k_z为Z轴磁强计标 度因数，b_xo为卫星轨道坐标系X轴磁场强度，b_zo为卫星轨道坐标系Z轴磁场强 度,k_x、k_z均为常数；

3)、采集N组测量数据组成测量模型：

$\left(\begin{array}{c}{H}_1\\ ·\\ ·\\ ·\\ H_N\end{array}\right)E=\left(\begin{array}{c}{G}_1\\ ·\\ ·\\ ·\\ G_N\end{array}\right)---\left(15\right)$

4)、采用最小二乘估计法估计得出磁强计零偏误差向量E：

E＝(H^TH)^-1H^TG      (17)

式(17)中，H为量测矩阵，G为量测矢量，H^T为H的转置；

5)、定义函数f：

$\left(\begin{array}{c}f(k_x,k_z,b_x^{ϵ},b_z^{ϵ})=\\ k_z{\tilde{b}}_x(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )-k_z{b}_x^{ϵ}(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )+\\ k_x{\tilde{b}}_z(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )-k_x{b}_z^{ϵ}(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )\end{array}\right)---\left(18\right)$

6)、以步骤2)中的标度因数和步骤4)得到的零偏误差为起始点， k_x,k_z,小量变化，采用N组测量数据对函数f进行计算，将每组等式记为 f_i，即得到：

F＝f₁²+f₂²+…+f_N²      (19)

式(19)中,F为N组量测数据构成的代价函数；

7)、以F作为标定的目标函数，并将优化算法设置为最小值搜索，在设定 标度因数和零偏参数范围内，令F为最小值的误差参数需要求解的量。

所述步骤1)中卫星俯仰角θ采用k时刻俯仰角θ(k)，其具体步骤为：

11)、利用太阳敏感器计算卫星的俯仰角：

${\theta }_s=\arctan \left\{\frac{-r_x{r}_{\mathrm{zo}}+r_z{r}_{\mathrm{xo}}}{r_x{r}_{\mathrm{xo}}+r_z{r}_{\mathrm{zo}}}\right\}---\left(10\right)$

式(10)中，θ_s为利用太阳敏感器输出计算得到的卫星俯仰角，r_x为测量的 卫星本体系X轴向太阳分量，r_z为测量的卫星本体系Z轴向太阳分量，r_xo为轨道 参考系下的X轴向太阳分量，r_zo为轨道参考系下的Z轴向太阳分量；

12)、得到Y方向陀螺输出与卫星的俯仰角速度的关系：

${\stackrel{·}{\theta }}_y={\omega }_y-ϵ---\left(11\right)$

式(11)中，为陀螺测得的卫星俯仰角速度，ω_y为陀螺输出，ε为陀螺输 出误差；

13)、综合式(10)与(11)，得出卫星的k时刻俯仰角θ(k)：

$\theta \left(k\right)=\frac{1}{2}{\theta }_s\left(k\right)+\frac{1}{2}(\theta (k-1)+\Delta ·{\stackrel{·}{\theta }}_y\left(k\right))---\left(12\right)$

式(12)中，θ_s(k)为k时刻太阳敏感器计算出的俯仰角信息，为k时刻 俯仰角速度信息，Δ为k-1到k时刻采样步长。

本发明的有益效果在于：(1)、本发明通过同时计算磁强计的零偏和标度因 数误差，并在姿态解算中予以补偿，从而有效地提高卫星导航精度；(2)、本发 明不需要专门针对标定算法设计新的滤波器，也不需要实时计算，有效地降低了 卫星运行成本和负载。

附图说明

图1为本发明标定流程示意图；

图2为俯仰角速率曲线图；

图3为俯仰角曲线图；

图4为俯仰角放大曲线图。

具体实施方式

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

如图1所示，本发明首先在光照区利用太阳敏感器和陀螺计算得到卫星的俯 仰角，然后以地面标定或上一次测量的标度因数误差为标准，只标定零偏误差， 最后重新建立多误差标定模型，采用最佳拟合的方法重新标定出零偏和标度因数 误差。在光照区以太阳敏感器与陀螺输出作为俯仰滤波器观测量，估计出卫星俯 仰角度和角速度，利用滤波输出量对磁强计误差进行估计，估计的结果进入滤波 器对磁场输出进行测量修正。

一、卫星姿态估计模型分析

1、磁强计误差模型

卫星磁强计的误差主要考虑为零偏和标度因数误差。剩磁是星上磁强计的主 要误差来源，剩磁造成的磁强计误差虽然是时变的，但是短期内可以近似认为是 常值，反映在误差项上即为星上磁强计零偏误差，同时考虑磁强计标度因数，则 误差模型如下：

$\tilde{b}=k·b+b^{ϵ}+\sigma ---\left(1\right)$

式(1)中，b为磁强计真值，k为磁强计标度因数，为磁强计测量值，b^ε为磁强计零偏误差，σ为磁强计随机量测噪声；

2、采用磁强计输出信息为量测量的俯仰滤波器

卫星的测量信息估计方法采用俯仰滤波器。卫星在俯仰Y轴方向上的姿态动 力学与X、Z轴解耦，则有：

$\stackrel{··}{\theta }=N_t/I_{\mathrm{yy}}---\left(2\right)$

式(2)中，θ是卫星相对于轨道的Y向姿态角，Y向一般被称为俯仰方向， 因此该姿态角简称卫星俯仰角，N_y为作用在Y轴上的控制力矩，I_yy为Y轴方向 的转动惯量；

状态方程：

$\left(\begin{array}{c}\theta (k+1)\\ \stackrel{·}{\theta }(k+1)\end{array}\right)=\left(\begin{array}{cc}1& \mathrm{\Delta t}\\ 0& 1\end{array}\right)\left(\begin{array}{c}\theta \left(k\right)\\ \stackrel{·}{\theta }\left(k\right)\end{array}\right)+\left(\begin{array}{c}\frac{{\mathrm{\Delta t}}^2}{2I_{\mathrm{yy}}}{N}_y\left(k\right)\\ \frac{\mathrm{\Delta t}}{I_{\mathrm{yy}}}{N}_y\left(k\right)\end{array}\right)---\left(3\right)$

测量方程：

$\theta =\arctan \left\{\frac{-b_x{b}_{\mathrm{zo}}+b_z{b}_{\mathrm{xo}}}{b_x{b}_{\mathrm{xo}}+b_z{b}_{\mathrm{zo}}}\right\}---\left(4\right)$

式(3)、(4)中，Δt为星载计算机滤波步长，k表示时刻，b_x为X轴磁强 计真值，b_z为Z轴磁强计真值，b_xo为卫星轨道坐标系X轴磁场强度，b_zo为卫星 轨道坐标系Z轴磁场强度。

二、基于多参数的磁强计迭代误差估计

1、磁强计多参数误差标定分析

由式(4)可得：

b_x(b_xosinθ+b_zocosθ)+b_z(b_zosinθ-b_xocosθ)＝0      (5)

考虑磁强计测量误差主要为零偏误差和标度因数误差，：

$b_x=({\tilde{b}}_x-b_x^{ϵ})/k_x,b_z=({\tilde{b}}_z-b_z^{ϵ})/k_z---\left(6\right)$

式(5)、(6)中，b_x为X轴磁强计真值，为X轴磁强计测量值，由磁强计 测量输出获得，为X轴磁强计零偏误差，k_x为X轴磁强计标度因数，由地面标 定或者上一次测量的标度因数误差暂定为常数，b_z为Z轴磁强计真值，为Z轴 磁强计测量值，由磁强计测量输出获得，为Z轴磁强计零偏误差，k_z为Z轴磁 强计标度因数，由地面标定或者上一次测量的标度因数误差暂定为常数，b_xo为 卫星轨道坐标系X轴磁场强度，由卫星预存的磁场表查表获得，b_zo为卫星轨道 坐标系Z轴磁场强度，由卫星预存的磁场表查表获得,所以有：

$\left(\begin{array}{c}({\tilde{b}}_x-b_x^{ϵ})(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )/k_x+\\ ({\tilde{b}}_z-b_z^{ϵ})(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )/k_z=0\end{array}\right)---\left(7\right)$

对式(7)进一步推导：

$\left(\begin{array}{c}{k}_z{\tilde{b}}_x(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )-k_z{b}_x^{ϵ}(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )+\\ k_x{\tilde{b}}_z(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )-k_x{b}_z^{ϵ}(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )=0\end{array}\right)---\left(8\right)$

写成矩阵形式：

$[{\tilde{b}}_x(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )-(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta ){\tilde{b}}_z(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )-(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )]\left(\begin{array}{c}{k}_z\\ k_z{b}_x^{ϵ}\\ k_x\\ k_x{b}_z^{ϵ}\end{array}\right)=0---\left(9\right)$

该方程为多元齐次方程，存在无穷多个解，因此不能直接用最小二乘的算法 来采用多组测量值求解。

2、磁强计多参数人迭代误差标定

由于磁强计的输出耦合了误差，俯仰滤波器中的姿态估算存在较大误差。在 光照区间太阳敏感器可用，陀螺也可短期开机，因此可以选择精度较高的太阳敏 感器和陀螺估计出姿态信息。

利用太阳敏感器信息计算卫星的俯仰角为：

${\theta }_s=\arctan \left\{\frac{-r_x{r}_{\mathrm{zo}}+r_z{r}_{\mathrm{xo}}}{r_x{r}_{\mathrm{xo}}+r_z{r}_{\mathrm{zo}}}\right\}---\left(10\right)$

式(10)中，θ_s为利用太阳敏感器输出计算得到的卫星俯仰角，r_x为测量的 卫星本体系X轴向太阳分量，r_z为测量的卫星本体系Z轴向太阳分量，r_xo为轨道 参考系下的X轴向太阳分量，r_zo为轨道参考系下的Z轴向太阳分量；

Y方向陀螺输出与卫星的俯仰角速度的关系为：

${\stackrel{·}{\theta }}_y={\omega }_y-ϵ---\left(11\right)$

式(11)中，为陀螺测得的卫星俯仰角速度，ω_y为陀螺输出，ε为陀螺输 出误差(含轨道角速率)。

综合式(10)与(11)，忽略误差项的影响，得到k时刻卫星俯仰角估计值为：

$\theta \left(k\right)=\frac{1}{2}{\theta }_s\left(k\right)+\frac{1}{2}(\theta (k-1)+\Delta ·{\stackrel{·}{\theta }}_y\left(k\right))---\left(12\right)$

式(12)中，θ(k)为k时刻俯仰角估计值，θ_s(k)为k时刻太阳敏感器计算出 的俯仰角信息，为k时刻俯仰角速度信息，Δ为k-1到k时刻采样步长。

分析磁强计的零偏和标定因数的特点，标度因数是磁强计测量与真实输出的 比例系数，变化量级相对较小，而零偏受星上剩磁影响严重，变化量级较大。当 考虑标度因数k_x,k_z为常数时，将方程(7)改写成以下零偏误差粗标定模型：

$\left(\begin{array}{c}\left[k_z(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )k_x(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )\right]\left(\begin{array}{c}{b}_x^{ϵ}\\ b_z^{ϵ}\end{array}\right)\\ k_z{\tilde{b}}_x(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )+k_x{\tilde{b}}_z(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )\end{array}\right)---\left(13\right)$

即：

H_iE＝G_i      (14)

式(14)中，H_i为量测量和已知量计算得到的用于标定的矢量， H_i＝[k_z(b_xosinθ+b_zocosθ)k_x(b_zosinθ-b_xocosθ)]；E为磁强计零偏误 差向量，G_i为量测量和已知量计算得到的用于标定的标量， $G_i=k_z{\tilde{b}}_x(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )+k_x{\tilde{b}}_z(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta );$

采集N组测量数据，测量值组成测量模型：

$\left(\begin{array}{c}{H}_1\\ ·\\ ·\\ ·\\ H_N\end{array}\right)E=\left(\begin{array}{c}{G}_1\\ ·\\ ·\\ ·\\ G_N\end{array}\right)---\left(15\right)$

则：

HE＝G      (16)

式(15)、(16)中，H为量测矩阵，G为量测矢量；

利用最小二乘估计法，可以估计出磁强计零偏误差向量E：

E＝(H^TH)^-1H^TG      (17)

式(17)中，H^T为H的转置；

将式(7)定义为函数f：

$\left(\begin{array}{c}f(k_x,k_z,b_x^{ϵ},b_z^{ϵ})=\\ k_z{\tilde{b}}_x(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )-k_z{b}_x^{ϵ}(b_{\mathrm{xo}}\sin \theta +b_{\mathrm{zo}}\cos \theta )+\\ k_x{\tilde{b}}_z(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )-k_x{b}_z^{ϵ}(b_{\mathrm{zo}}\sin \theta -b_{\mathrm{xo}}\cos \theta )\end{array}\right)---\left(18\right)$

以地面标定或上一次测量的标度因数和标定后的零偏误差为起始点， k_x,k_z,小量变化，采用多组测量数据对函数f进行计算；在多组测量数据 下，可得到不同的函数计算值，将每组等式记为f_i。

因此，记

F＝f₁²+f₂²+…+f_N²      (19)

式(19)中,F为N组量测数据构成的代价函数；

将F作为标定参数问题的目标函数，并将优化算法设置为最小值搜索。则在 设定标度因数和零偏参数范围内，其中标度因数范围可以取[标度因数起始点 -0.1，标度因数起始点+0.1]，零偏范围可以取[零偏起始点-0.1，零偏起始点 +0.1]，令F为最小值的误差参数为需要求解的量，即k_x,k_z,

在光照区，利用式(17)，其中俯仰角计算采用式(12)估计值，轨道系磁 场强度通过轨道信息查表求得，则可以估计出磁强计误差。由于此时估计的角度 和角速度信息与剩磁误差不相关，因此计算出的俯仰角和俯仰角速度不仅可以用 于姿态控制，也可以作为辅助信息进入滤波器对磁场输出进行测量修正。

3、仿真验证

卫星轨道为10∶30AM太阳同步轨道,高度496km。磁强计量测噪声均方 差30nT，输出频率10Hz。设置卫星的剩余磁矩为[0.30.30.3]A·m²，为 更进一步模拟磁场误差，假设磁强计标度因数为0.988，零偏误差为800nT，地 面标度因数试验结果为0.991。卫星转动惯量为diag[0.80.80.8]kg·m²，姿 态控制系统使卫星保持对地稳定，并加入了剩磁扰动力矩，其余干扰力矩假设为 白噪声，均方差为5×10^-7N·m。太阳敏感器视场范围±30°×±30°，测量精 度优于0.5°(3σ)，陀螺的测量精度约为0.003(°)/s。

每次进入光照区开始数据采集，同时进行实时估计，估计出的误差在滤波器 中用于量测修正。假设标度因数为0.991的常量，经过两个轨道周期后的数据采 集，X轴磁强计零偏误差结果为988nT，Z轴磁强计零偏误差估计结果为637T。 然后开始误差的迭代标定，标度因数以0.001为步长，迭代范围为[0.991-0.1, 0.991+0.1]，零偏以30nT为步长，迭代范围为分别为[988-600,988+600]， [637-600,637+600]。重新经过两个轨道周期的数据采集计算后，k_x,k_z,标 定结果分别为0.986、0.990、904、723。

如图2、3所示,将只标定磁场零偏与同时标定标度因数和零偏两种情况下的 参数用于滤波补偿，补偿后的俯仰角和俯仰角速度曲线，可见在两种补偿方式下 系统都是稳定的。

如图4所示,将俯仰角曲线局部放大，b曲线为只标定磁场零偏的补偿结果， a曲线为同时标定标度因数和零偏的补偿结果,可以看到a曲线相比b曲线有明 显改善，整体波动更为平缓。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技 术人员来说，在不脱离本发明原理的前提下还可以做出若干改进，这些改进也应 视为本发明的保护范围。