在过去的15年内，随着不同应用领域任务需求的日益增加，小型化、高性价比卫星技术得到了蓬勃发展。小型卫星能借助电子产业大规模工业生产的基础和架构，以尽可能小的质量、体积和批量化生产的方式实现任务需求，主要面向教学与科研、低轨通信、新技术验证，以及未来空间遥感组网、空间碎片监测等任务^[1]。目前国际认同的对小型卫星依据质量的分类如下：微小卫星, 质量为10~100 kg; 纳卫星, 质量为1~10 kg; 皮卫星, 质量为0.1~1 kg; 飞卫星, 质量为0.01~0.1 kg。小型卫星正朝着更小、更轻、更廉价的方向发展^[2]。文献[2]指出和1U的立方星相比，飞卫星具有成本低、通过冗余缓解固有风险以及组网控制等单个纳卫星难以具备的优势。文献[3]同样表明相对单个大卫星而言，使用飞卫星组建的分布式航天器系统具有更好的灵活性和鲁棒性。

目前对于飞卫星的任务设计^[4]、系统设计^[2]、能耗设计^[5]和传感器设计^[6]等领域的研究已经逐渐丰富。文献[7]提出了一种使用MEMS控制力矩陀螺(control moment gyroscope, CMG)设计的飞卫星姿态控制执行器, 对MEMS CMG进行建模并详细介绍了一个的飞卫星设计方案。为了降低制造成本, 文献[2]采用商业化元件选型, 对嵌入式卫星的质量、能源、硬件、软件、散热问题做了详细的介绍, 但是对姿态控制系统分析不足, 未能针对芯片的性能提出软件设计方案。文献[8]重点分析了PCBSat和WikiSat 2颗飞卫星在功能、结构、功耗等方面的特点。PCBSat设计有2个太阳敏感器、1个GPS导航接收机和被动气动控制系统, 而WikiSat包含4个光学传感器、1个三轴加速度计、1个三轴陀螺仪以及磁力距器用于姿态控制。但是由于飞卫星的质量、体积和功耗方面的限制, 目前飞卫星很少设计姿态和轨道控制任务^[4]。

与纳卫星和皮卫星姿态控制相比, 飞卫星姿态控制系统的MCU计算能力偏弱, 且存储空间更小。例如，现有卫星PCBSAT^[8]采用了ATmega 128L微控制器，FLASH容量为128 kB；RyeFemSat^[9]飞卫星采用了32 kB FLASH的CC2510控制器；WikiSat^[8]飞卫星所采用的控制器ATmega328P也仅有32 kB FLASH。因此考虑工程实现的可行性, 飞卫星姿态控制建模也应该与纳卫星和皮卫星姿态控制建模不尽相同^[7]。工程上一般使用历表的方法解决此类复杂模型计算问题^[10]。本文根据飞卫星处理器计算性能较弱的特性, 对卫星姿态控制系统中计算复杂度较高、存储空间较大的地磁场模型进行简化, 设计了均匀模型和非均匀模型, 并通过将模型应用到卫星姿态确定中对比分析, 最后总结并选取适用于飞卫星MCU的80点非均匀地磁场历表模型。

1 均匀地磁场历表模型 1.1 地磁场标准模型

通常可将地磁模型分为高空全球模型和区域模型, 全球模型主要以国际地磁参考IGRF为代表, 其模型的建立利用地面、海洋、高空和卫星地磁测量数据, 通过Gauss理论而建立。自1968年起, 国际地磁与高空物理协会(IAGA)每隔5年会发布国际地磁参考模型(IGRF)。

本文的地磁场模型采用IAGA提供的IGRF2010模型。使用球谐波函数表示的地球磁位势为

(1)

式中：R_e为地球半径; r为格林尼治起算的东经; λ和θ分别为地理经度和地心余纬; g_n^m和h_n^m为基本磁场的高斯系数; P_n^m(cosθ)为n次m阶缔合勒让德函数。

地球磁势位在三轴方向的梯度即为地磁场矢量:

(2)

式中，P_n^m(θ)和使用递归方式计算。

仿真分析所用到的地磁矢量由IGRF模型模拟得到, 但是受限于星载计算机计算速度和星上存储容量, 实际星上加载地磁模型的级数不能取过高^[11]。为了本文的仿真分析, 本文选择了10阶球谐级数的地磁场模型, 模拟了轨道高度为650 km, 轨道倾角为97°, 升交点赤经10°的卫星轨道上的地磁场矢量随时间变化曲线如图 1所示。

1.2 地磁场历表模型

地磁场模型只有在精度比较高时, 才能获得满意的精度, 这种实时迭代解算的计算量是星上计算机很难承受的。以IGRF2010模型为例, 该模型包含2003年最新修订的195个高斯系数, 在弱太阳活动时间段内, 计算精度可达100 nT以下^[12]。历表模型以损失精度、增加存储量为代价, 保证了地磁场模型的计算效率和双矢量算法的正常运行。

本文设计地磁场历表模型的主要思想为:在存在有效GPS数据的情况下, 通过卫星的轨道位置(GPS)推算轨道6根数中的真近点角, 进而利用真近点角和地磁场强度的对应关系进行建表和查表, 然后对查到的历表数据进行线性插值, 最终解算出卫星所处轨道位置的磁场强度数据。算法设计流程如图 2所示。在没有有效GPS数据的情况下, 通过卫星运行轨道推算卫星轨道位置的真近点角, 然后查表、插值得到卫星所处轨道位置的磁场强度。

1.3 历表设计

本文设计的地磁场历表模型充分考虑了地磁场模型在不同应用场景下可能存在的不同精度需求, 即考虑了同一轨道不同表格长度的历表模型。本文设计的地磁场历表模型的存储格式由三部分组成:起始真近点角θ₀、表格长度L和地磁场数据M₀。图 3表述了不同起始真近点角和表格长度的地磁场历表模型下的表格存储格式。

本文存储实数数据采用单精度浮点格式, 因此存储初始真近点角和表格长度均使用4字节, 而存储轨道上每一点的三轴地磁场需要占用12字节的存储空间, 因此存储L表格长度为的历表模型需要12L+8字节。

1.4 选取表格长度

卫星轨道上的地磁场是连续分布的, 理论上来说可以采集到无穷多点的地磁场数据, 但是实际上表格长度选取的过大无疑会增加存储负担, 同时冗余数据的存储也无助于提高算法精度。相反, 表格数据过小的情况下, 查表法历表模型本身存在较大的误差甚至错误, 导致无法利用这种历表模型对姿态确定的双矢量算法精度进行测试。因此选择一个合适的表格长度显得尤为重要。

1.5 解算真近点角

在卫星飞行的过程中, 认为卫星的运行轨迹为圆形, 轨道高度固定, 其圆点在地心。轨道6根数是经典万有引力定律描述天体按圆锥曲线运动时所必需的6个参数, 如图 5所示。

轨道上任意一点的GPS数据均可由轨道6根数计算确定, 轨道6根数确定卫星经度和纬度的计算如公式(3)所示。

(3)

式中：λ₀为0时刻的升交点赤经；λ_s(t)为经度；φ_s(t)为纬度；i为轨道倾角；θ为真近点角；ω_e为地球自转角速度。

当卫星所处的轨道和轨道位置即λ_s(t)和φ_s(t)的值已知时, 可以通过(3)式求解出真近点角θ的值, 真近点角的取值范围是0~360°。

(4)

在未给出GPS数据的情况下, 卫星的运动可近似看成匀速圆周运动, 即卫星绕地心飞行的角速度ω_s为定值, 真近点角θ随时间均匀变化。因此, 在给定初始真近点角θ和飞卫星飞行时间Δt的情况下, 可以推算出当前时刻卫星所处轨道位置的真近点角为

(5)

使用推算真近点角的方法, 只需要卫星轨道参数、飞行角速度和至少一个初始真近点角或者GPS数据, 即可推算之后任意时刻卫星所处的轨道位置和该处地磁场强度的大小。

2 非均匀地磁场历表模型

观察图 1发现磁场数据并不是随着时间均匀变化, 而是在某些时间段磁场曲线弯曲程度大, 另外一些磁场数据弯曲程度小, 因此采用均匀采样的方法势必造成曲率小的时间段内磁场数据冗余采样, 而在曲率较大的时间段内磁场数据不足造成精度损失。为了解决磁场曲线弯曲程度和采样点数之间的矛盾, 本文借用曲率的概念, 在曲线曲率大的地方多采样, 曲率小的地方少采样。地磁场曲线的曲率计算过程中需要对公式(2)求二阶微分, 因此计算复杂度较高, 本文采用离散曲率的方法代替直接对地磁场直接求导的过程, 从而实现计算过程的简化。

2.1 极坐标归一化

考虑到磁场曲线的横坐标表示时间, 纵坐标表示磁场强度, 直接对图 1曲线求曲率是没有意义的, 因此在求曲率之前先把磁场曲线转换成极坐标形式。这种使用一维曲线来表示二维目标轮廓的方法在形状特征提取方面也有应用^[13]。

将图 1转换成极坐标的形式考虑到三轴磁场曲线是相互独立的, 于是再对每轴进行单独的归一化。需要注意的是, 必须将极坐标下归一化磁场曲线转换成直角坐标系下归一化地磁场曲线的形式之后, 才可以对曲线求导。其原因在于直角坐标转换成极坐标的过程中图形不同区域曲率的相对大小发生了改变, 损失了曲率的相对信息。

2.2 U弦长曲率

U弦长曲率是一种离散曲率计算方法^[14], 与现有的离散曲率计算方法相比, U弦长曲率具有更强的抗旋转性和抗噪性, 适用于完成曲线匹配等对曲率计算稳定性要求高的一类任务。方法基本思想是:对于输入的参数U, 按照欧氏距离在曲线当前点处确定一个支持领域, 并应用文献[15]中的曲线精化策略改进计算精度, 由此计算离散曲率。

记:l={p_i:(x_i, y_i, i=1, 2, …, N)}为一条数字曲线, 计算p_i处的U弦长曲率时, 应用与支持领域前后臂矢量夹角相关的一个余弦值作为离散曲率, 具体计算公式为

(6)

式中，Q_i是p_i点处的离散曲率; s_i=sign[(x_i-x_i^b)(y_i^f-y_i^b)-(x_i^f-x_i^b)(y_i-y_i^b]为离散曲率值的符号; D_i=p_i^fp_i^p是p_i^f、p_i^b这两点间的欧氏距离; (x_i^f, y_i^f), (x_i^b, y_i^b)分别是p_i^f和p_i^b的坐标; U为支持邻域的长度, 具体参见文献[11]。

2.3 K曲率加权非均匀采样

K曲率加权的方法是本文提出的一种对每个离散点进行加权的方法, K代表离散点的固定权重, 通过实验的方法选取最优的K值。令地磁场数据为1列N个有序的离散点, 在p_i处的离散曲率为Q_i, 则所有离散点曲率模的和为

(7)

对N个有序离散点的离散曲率进行归一化处理, 得到p_i点处的离散曲率模为|Q_i|/Q_sum。在对地磁场数据进行非均匀离散化处理时, 若只考虑每个点的曲率权重, 则容易造成所采样的点过度地聚集到曲率较大的位置处。所以对于每个离散点除了考虑曲率权重外, 还应附加自身固有的权重, 用以平衡这种过度聚集的情况。本文将固有权值K赋值给N个点, 则p_i点的权重为K/N。因此, K曲率加权的采样方法需要满足

(8)

式中, n为采样点数, j=1, 2, …, n, m_j表示采样点位置对应标准模型中的位置。

易知, K的选取决定了曲线的非均匀采样结果, 当K无限接近于正无穷时, K曲率加权的采样方法便蜕化成了均匀采样。K曲率加权采样方法是曲率和均匀采样的融合, 不同的K值会得到不同的非均匀采样结果, 通过实验方法筛选出最优的K, 使得非均匀采样的点能够更好地描述原曲线。

图 6是采用K曲率加权非均匀采样方法得到的80点地磁场模型, 和文章的设计意图相同, 采样密度较高的地方集中在曲率较大的区域。可以看出随着采样点数的减小, 三轴地磁场图形变得越来越简陋, 离散点所表达的信息也越来越少, 图形的轮廓变得越来越模糊。

在K曲率加权采样选取了采样点之后, 使用线性内插获得地磁场数据l_K={B_i^K:(1=1, 2, …, N)}, 与标准地磁场数据l_S={B_i^S:(i=1, 2, …, N)}之间的平均误差记作

(9)

本文选取的一个轨道的地磁场数据的N=5 875。E用于描述2条磁场曲线的偏离程度, E越小, 表示2条曲线越贴近。

图 7表示采样点数为80时K和E的关系, 从图中可以看出, 采样点数相同时, 均匀采样的E为常数, 记作E^J, 而K曲率加权采样的E随K变化, 取最小的E记作E_min^K。总能找到存在E_min^K < E^J, 表示K曲率加权采样得到的磁场曲线和原磁场曲线更加贴近, 从而证明了算法的有效性。另外, 随着采样点数的减小, E_min^K和E^J都在增加, 表示采样点数减小, 2种方法的拟合曲线和原曲线之间的误差都增加。

3 飞卫星姿态确定实验 3.1 双矢量姿态确定模型

本文中飞卫星的姿态确定系统采用太阳敏感器和磁强计的组合作为敏感器, 并基于矢量算法进行飞卫星姿态确定算法的设计。双矢量定姿的思路如下:

1) 根据轨道坐标系中的单位矢量S₀, B₀, 在轨道坐标系中建立新的正交坐标系L, 各坐标表轴的单位矢量为:

(10)

2) 根据卫星本体坐标系中的单位矢量S_b, B_b在卫星本体系中建立一个正交坐标系S, 各坐标轴的单位矢量为:

(11)

3)定义矩阵M_L=[l₁   l₂   l₃], M_S=[s₁    s₂   s₃], 可得M_s=A_b0M_L, 则存在唯一的正交姿态矩阵, 满足

(12)

图 8描述了双矢量姿态确定的流程, 其中对地磁场数据的使用包括图中虚线框内的3个流程。首先, 测定轨道位置、加载轨道系地磁场数据; 其次, 初始化卫星的轨道位置并进行轨道推演; 最后使用地磁场模型解算对应的轨道系地磁场矢量数据。图中磁场强度模型和太阳能电池板电压模型引自文献[16]。

3.2 蒙特卡洛模拟

在模拟地磁场数据和太阳矢量数据的过程中, 由于随机误差的引入, 往往导致模拟的量测数据和实际量测的数据之间存在偏差, 为了降低此类偏差以提高模拟精度, 本文采用蒙特卡洛模拟降低此类误差的影响。蒙特卡洛方式是一种概率统计理论为指导的数值计算方法, 指的是使用随机数来解决很多计算问题的方法。

本文使用蒙特卡洛模拟执行M次单轨姿态确定解算, 将IGRF地磁场模型作为轨道系地磁场数据解算出的欧拉角记作Q_{i, j}^S, 其中i=1, 2, …, N, j=1, 2, …, M。同样, 使用均匀采样地磁场模型解算的欧拉角记作O_{i, j}^U, 使用K曲率加权地磁场模型解算的欧拉角记作O_{i, j}^K。为描述O_{i, j}^U和O_{i, j}^K分别与O_{i, j}^S之间的误差, 本文引入了平均均方根误差θ_ARMSE的概念, 如公式(13)~(14)所示。其中, 上标SU表示均匀采样, 上标SK表示非均匀采样。θ_ARMSE越小说明O_{i, j}^U或O_{i, j}^K分别与O_{i, j}^S之间的差异越小, 即采样结果越精确。

(13)

(14)

3.3 实验结果和分析

实验中, 选取蒙特卡洛模拟次数M=50, 轨道周期N=5 875 s, 地球半径为6 385 km, 轨道高度为650 km, 升交点赤经为10°, 近地点幅角为10°, 离心率为0.1, 轨道倾角为97°, 轨道常量为3.986×10¹⁴, 轨道角速度由轨道常量和轨道半径计算得出; 地磁场模型噪声的均值为0, 标准差为10^-8 T; 太阳矢量模型的均值为0, 标准差为0.01。

这里以X轴为例, 从图 9中可以看出θ_ARMSE^SK < θ_ARMSE^SU始终成立, 说明使用K曲率加权的地磁场模型得到的姿态确定结果和使用标准姿态确定结果更加接近, 从而证明了K曲率加权的方法比均匀采样的方法更加优越; 另外, 当采样点数较少时, θ_ARMSE^SK和θ_ARMSE^SU之间差异更明显, 说明采样点数越少, K曲率加权采样的优点突出的更明显; 最后, 当采样点数为80左右时, 三轴θ_ARMSE^SK均趋于稳定, 因此在不影响姿态确定精度的前提下, 本文选取80点作为地磁历表模型的采样个数较为合理。

图 10为使用80采样点的K曲率加权地磁场模型得到的姿态确定结果，从图中可以看出其姿态确定精度维持在7°以内。

3.4 应用前景分析

本文的研究表明，对于一条轨道(5 875 s)三轴地磁场数据，使用80采样点的地磁场历表模型基本可以替代IGRF模型在姿态确定中的作用。文献[10]提出了网格化地磁场模型，即使用经纬度数据去映射存储磁场数据，若以4字节的float格式存储磁场数据，则共需要128.2 kB。由于飞卫星没有变轨功能，因此网格化磁场模型中存储了大量的数据冗余。相比之下，若IGRF模型和80采样点的地磁场数据以4字节的float格式存储，则各需要70.5 kB和0.96 kB。由此可见，本文提出的基于K曲率加权的地磁场历表模型不仅占用存储空间较小，满足在前文提及的3种飞卫星的FLASH中轻量化存储的需求，同时可实现将FLASH中的地磁场数据一次性加载到MCU的RAM中，从而达到对地磁数据快速访问的目的。所以本文设计的磁场历表模型不仅极大程度上节约了存储空间，而且提高了程序的运行速度，从而缩短了算法的控制周期。

4 结论

本文致力于解决地磁场模型简化计算量和存储空间的问题，给出了地磁场历表模型的表格存储格式；并提出了K曲率加权的非均匀采样地磁场模型，在采样点数固定的情况下，给出非均匀点的选取方法；最后，为验证K曲率加权模型的正确性，文章采用蒙特卡洛模拟方法对卫星姿态确定过程进行分析和验证。K曲率加权模型在满足卫星姿态确定精度的前提下，尽量减少选取的采样点数，从而降低所需的存储空间，在计算能力偏弱、存储空间较小的飞卫星姿态确定领域具有非常旷阔的应用前景。