本发明属于航天器控制领域，涉及一种航天器姿轨控系统中星敏感器系统误差的在轨辨识与实时补偿方法。
背景技术：
星敏感器的工作原理是通过对星空进行拍照，提取星点与其内置的星库进行比较，获得星敏感器在惯性空间中的高精度姿态信息。星敏感器本身的测量精度很高，作为航天器姿轨控系统的关键测量部件，其测量精度直接决定了航天器平台的姿态确定精度。
但是，星敏感器属于舱外单机，其光照和热控环境比较恶劣，受制于星敏感器安装面的热变形、星敏感器本身结构的热变形以及星图匹配因素等影响，星敏感器测量的系统误差主要表现为轨道周期不规则交变的形式。由于航天器轨道在惯性空间中缓慢进动，航天器光照条件、航天器和星敏感器的热变形、星敏感器视场扫描星空区域以及轨道周期本身等均是缓慢变化的，因此，星敏感器的轨道周期系统误差也是缓慢变化的。要提升航天器姿态测量精度，需对星敏感器轨道周期系统误差进行精确的补偿。
陀螺测量精度非常高，一般陀螺短时间内积分误差远高于星敏测量精度。但由于工艺等因素影响，陀螺不可避免存在常值漂移系统误差，这大大约束了采用高精度陀螺积分作为航天器姿态基准的应用前提。
采用高精度陀螺作为航天器基准的前提是将陀螺常值漂移系统误差控制在一个极小的范围内。虽然星敏感器存在轨道周期系统误差，但在相同的轨道位置，星敏感器测量系统的误差变化极小(可近似认为是0)，基于此原理，可以消除星敏感器修正陀螺时由于星敏感器热变形等系统误差带来的陀螺常值漂移修正误差。
此外，陀螺的安装约束比星敏感器少很多，对于某些精度要求较高的载荷，可将高精度陀螺与平台高精度载荷共基准安装，将由于热控和热变形等因素对姿态指向精度的影响减至最小，因此高精度陀螺更适合作为航天器的姿态基准。
目前还没有采用高精度陀螺对星敏轨道周期系统误差进行在线补偿的公开技术。现有公开文献中，专利《一种基于星敏的地平仪系统误差在线补偿的方法》(授权号：ZL201310266392.1)中，以星敏感器为基准，系统可以在线识别地平仪轨道周期系统误差并在线补偿，但要对星敏感器轨道周期系统误差进行补偿，需要精度更高的姿态基准，此专利中并未涉及。现有对星敏感器轨道周期系统误差补偿，普遍的做法是将陀螺和(或)载荷数据下传至地面，地面经过事后处理识别出星敏感器的系统误差，并将此系统误差进行曲线拟合，再上注至星载计算机。这种做法弊端如下:1.不能在线实时识别并补偿星敏的轨道周期系统误差，对星敏感器轨道周期系统误差的精度低；2.需要地面干预，地面工作量大，可靠性相对较差。
技术实现要素：
本发明解决的技术问题是：克服现有技术的不足，提供了一种利用陀螺修正星敏感器轨道周期系统误差的方法，解决了星敏感器系统误差无法被在轨实时辨识及修正的技术问题，提升了星敏感器系统误差修正的实时性和精确性，减轻了地面人员的工作量。
本发明的技术解决方案是：一种利用陀螺修正星敏感器轨道周期系统误差的方法，包括如下步骤：
(1)构建补偿差值表；所述的插值补偿表用于存储星敏感器三轴姿态系统误差和航天器所处轨道相位的关系，所述的轨道相位采用轨道钟描述，轨道钟是控制系统星载软件步长ts的整数倍，轨道钟以轨道上的某个点为起点，以ts为单位进行累加积分，待积分到刚好一个轨道周期时轨道钟清零，重新开始累加；
(2)利用星敏感器的姿态对陀螺的常值漂移进行修正；
(3)采用陀螺积分值来辨识星敏感器轨道周期系统误差，并对补偿差值表进行更新；
(4)查找更新后的补偿差值表，将星敏感器三轴姿态直接扣除补偿差值表中对应的值，得到不含系统误差的星敏感器三轴姿态。
所述的轨道钟的计时方法为：
(11)轨道钟累加，torbit(k)＝torbit(k-1)+ts，其中torbit(k)为当前节拍轨道钟，torbit(k-1)为上节拍轨道钟；
(12)记录当前时刻航天器轨道幅角与180°之差Δu(k)＝u-180°，其中u为当前时刻航天器轨道幅角，航天器轨道幅角为航天器轨道六要素中的近地点幅角与真近点角之和；
(13)判断Δu(k)·ΔU<0是否成立，其中ΔU是记录的上一节拍航天器轨道幅角与180°之差，初始值为0；如果成立则轨道钟清零并进入下一步，如果不成立，则直接转至下一步；
(14)将当前时刻航天器轨道幅角与180°之差Δu(k)赋值给ΔU。
所述步骤(2)中的修正方法为：
(21)根据积分时间t(k)判断当前是否是陀螺积分第一节拍，如果是则跳至步骤(22)，如果否则跳至步骤(23)；
(22)记录当前星敏感器的姿态，记为θstr0并结束；
(23)判断积分时间t(k)是否到达轨道周期T，如果是则跳至步骤(25)，如果否则跳至步骤(24)；
(24)陀螺一个节拍积分θgyro(k)＝θgyro(k-1)+ω，其中θgyro(k)为当前节拍陀螺积分姿态，θgyro(k-1)是上一节拍陀螺积分姿态，ω是当前节拍陀螺测量角速度；时间积分t(k)＝t(k-1)+ts，t(k)为当前节拍陀螺积分时长，t(k-1)是上一节拍陀螺积分时长，然后跳至步骤(26)；
(25)记录当前星敏感器的姿态，记为θstrf；计算陀螺的常漂残量
计算结束跳至步骤(27)；
(26)判断陀螺积分是否超出阈值，如果是则跳至步骤(28)，如果否则结束；
(27)判断陀螺常漂残量Bremain是否小于阈值BC，如果是则跳至步骤(29)，如果否则跳至步骤(28)；
(28)清除陀螺修星敏标志flag＝0，然后跳至步骤(30)；flag是表示是否可以采用陀螺积分值来辨识星敏感器轨道周期系统误差的标志；
(29)建立陀螺修星敏标志flag＝1，然后跳至步骤(30)；
(30)计算陀螺常值漂移为本轨道周期计算得到的陀螺常值漂移，其中ktl∈[0,1]为陀螺常值漂移更新权系数，为上次计算得到的陀螺常值漂移；同时
积分时间清零t(k)＝0；
陀螺积分清零θgyro(k)＝0；
记录当前轨道钟T0＝torbit(k)；
记录陀螺积分残差θ0＝θgyro(k)-θstrf+θstr0。
所述步骤(3)中更新补偿差值表的方法为：
(41)判断flag＝1是否成立，如果成立则跳至步骤(42)；如果否则退出；
(42)判断当前轨道钟torbit(k)是否大于记录时刻轨道钟T0，如果torbit(k)-T0≥0，则跳至步骤(43)，如果否则跳至步骤(44)；
(43)修正量偏差计算计算结束跳至步骤(45)；Δθ为当前节拍基于陀螺积分修正星敏感器时，应扣除的修正量偏差；
(44)修正量偏差计算计算结束跳至步骤(45)；
(45)计算即将要更新的补偿表行号ΔT为ts的整数倍；
(46)判断当前节拍是否需要对相应的补偿表进行更新，如果满足torbit＝(M-1)·ΔT，则转至步骤(47)，否则退出；
(47)计算第M个点的平均修正量θxz＝θgyro(k)-Δθ-θstr(k)，对第M行数据进行更新其中θxz是当前节拍计算得到的星敏感器系统误差值，上方横杠表示取平均值，θgyro(k)是当前节拍陀螺积分值，θstr(k)为系统误差经补偿后的星敏感器实测值，kxm∈[0,1]，计算结束跳至步骤(48)；
(48)判断M＝1是否成立，如果成立则跳至步骤(49)，如果否则退出；
(49)将第一个点修正值同时存入最后一个点，计算结束退出。
所述步骤(4)中查找更新后的补偿差值表时，如果能够直接找到对应的轨道钟，则直接读取该轨道钟对应的误差值参与计算；如果无法直接查找到对应的轨道钟，则采用插值的方法获得参与计算的误差值，具体为：根据当前轨道钟torbit，当前轨道钟前一行补偿插值表中的内容Block[M]_fy，当前轨道钟后一行补偿插值表中的内容Block[M+1]，计算得到当前轨道钟对应的星敏感器系统误差：
本发明与现有技术相比的优点在于：
第一：本发明提出了一种基于陀螺的补偿星敏轨道周期系统误差的方法，本方法不需要载荷信息，仅在姿轨控内部即可完成星敏轨道周期系统误差的补偿，具有普适性；充分发挥星敏感器和陀螺的特点及优点，先采用在同一轨位的星敏感器姿态与陀螺积分一个轨道周期的姿态进行对比，消除由于星敏感器热变形带来的陀螺常值漂移修正误差，使陀螺常值漂移得到精确修正，最终使得陀螺积分精度满足辨识星敏感器轨道周期系统误差标定的要求；
第二：采用基于陀螺的星敏轨道周期系统误差在线辨识及自主补偿，可增加星敏轨道周期系统误差补偿精度。基于上述技术方法，采用星敏感器轨道周期系统误差在线辨识及自主补偿技术，辨识得到的星敏感器系统误差会被立即应用，具有非常好的时效性，有助于提高星敏感器轨道周期系统误差的提高；
第三：采用基于陀螺的星敏轨道周期系统误差在线辨识及自主补偿，可大幅减小地面工作量，增加星敏轨道周期系统误差补偿精度和补偿可靠性。同时，采用星敏感器轨道周期系统误差实时辨识更新技术，可以省去地面对在轨数据的处理，大幅减小地面工作量，增加星敏轨道周期系统误差补偿精度和补偿可靠性。
附图说明
图1为本发明的轨道钟工作流程图；
图2为本发明的陀螺常值漂移在线修正流程图；
图3为本发明采用陀螺修正星敏感器轨道周期系统误差的流程图；
图4为本发明采用插值算法实现星敏感器轨道周期系统误差补偿的流程图。
具体实施方式
第一部分：制定插值补偿表
插值补偿表用于存储星敏感器三轴姿态系统误差，为了减小星上数据存储量，每隔一段时间ΔT提取星敏感器三轴姿态系统误差，并存储在表的对应位置。在存储表对应的轨道相位点，将星敏感器三轴姿态直接扣除存储表的值即得到不含系统误差的星敏感器三轴姿态。在存储表以外的轨道相位点，先采用两点间线性插值算法(具体算法后面详述)计算得星敏感器当前轨道相位点的系统误差，再将星敏感器三轴姿态扣除此系统误差，得到不含系统误差的三轴姿态。
插值补偿表分为四列，第一列为轨道钟，体现了星敏感器所处的轨道相位点，第二列为星敏感器滚动姿态系统误差，第三列为星敏感器俯仰姿态系统误差，第四列为星敏感器偏航姿态系统误差。
轨道钟定义：以某个点(如降交点)为起点，以控制系统星载软件步长ts为单位进行累加积分，待积分到刚好一个轨道周期左右，轨道钟清零，重新开始累加。不难看出：轨道钟是ts的整数倍，轨道钟与星敏感器所处的轨道相位点具有一一对应的关系。
表格的行数主要取决于星敏感器三轴姿态系统误差特性，和系统对星敏感器三轴姿态系统误差补偿精度要求，一般来说，表格行数越多，就能更加真实还原星敏感器三轴姿态系统误差，但是占用星上空间资源也相应越大，具体表格行数根据实际情况定。每ΔT(ΔT为ts的整数倍)选取一个点，则最终插值补偿表的行数为：其中，floor()为截去小数取整函数。表格第一列轨道钟为ΔT的整数倍。
插值补偿表的内容为星敏感器三轴姿态系统偏差，其计算方法如下：
当轨道钟为ΔT的整数倍时，记录下星敏感器三轴姿态与陀螺积分三轴姿态之差，存入表格对应位置。
由于航天器轨道周期不是定值，如某航天器轨道周期T变化区间为[Tlow,Thigh]，因此星上无法根据航天器轨道钟是否到轨道周期来决定是否清零。相应地，轨道钟的定义是以航天器降交点为起点(轨道钟为0)，以控制系统步长ts为单位进行累加积分的，因此，轨道钟清零判断准则为：航天器轨道幅角第一次到达180°(降交点)。
如图1所示，轨道钟工作步骤如下：
步骤1：
轨道钟累加torbit(k)＝torbit(k-1)+ts，其中torbit(k)为当前节拍轨道钟，torbit(k-1)为上节拍轨道钟，ts为控制系统星载软件步长，计算结束转至步骤2；角标为k表示对应当前节拍，角标为k-1表示对应上一节拍。
步骤2：记录当前时刻航天器轨道幅角与180°之差：Δu(k)＝u-180°，其中u为当前时刻航天器轨道幅角，其定义为航天器轨道六要素中的近地点幅角与真近点角之和，属于领域内通用定义。计算结束转至步骤3；
步骤3：判断轨道钟是否满足清零条件，判据为：Δu(k)·ΔU<0？，其中ΔU是记录的上一节拍航天器轨道幅角与180°之差，软件默认状态为0，如果是，则转至步骤4(轨道钟清零)，如果否，则转至步骤5；
步骤4：轨道钟清零，计算结束，转至步骤5；
步骤5：将当前时刻航天器轨道幅角与180°之差Δu(k)赋值给ΔU，供下次判断轨道钟是否需要清零判断用。
因此，用于存储星敏感器三轴姿态系统误差的插值补偿表，是一个列数为4，行数为的表，其体现了星敏感器三轴姿态系统误差在一整轨的变化情况。
星敏感器三轴姿态系统误差为星敏姿态与真实姿态之差，由于无法得知真实姿态，但可以用高精度陀螺积分信息来代替真实姿态，这也是本发明的具体操作方式。
第二部分：星敏感器修正陀螺常值漂移算法：
本发明的核心思想，是以高精度陀螺积分为基准，来辨识及在线更新星敏感器三轴轨道周期系统误差，最终达到星敏感器轨道周期系统误差在线补偿的目的。要满足陀螺积分三轴姿态精度足够高的条件，必须对陀螺常值漂移进行精确修正。
在进行星敏感器修正陀螺常值漂移算法时，星敏感器三轴角度需满足小角度条件，此时航天器三轴姿态可视为相互解耦。本发明以其中一轴姿态θ为例进行设计说明，如图2所示：
步骤1：根据积分时间t(k)判断当前是否是陀螺积分第一节拍，如果是，则跳至步骤2；如果否，则跳至步骤3；
步骤2：记录当前星敏感器的姿态，记为θstr0，然后结束；
步骤3：判断积分时间t(k)是否到达轨道周期T(T是航天器轨道周期，一般根据航天器当前轨道计算得到)，如果是，则跳至步骤5；如果否，则跳至步骤4；
步骤4：陀螺一个节拍积分θgyro(k)＝θgyro(k-1)+ω，其中，θgyro(k)为当前节拍陀螺积分姿态，θgyro(k-1)是上一节拍陀螺积分姿态，ω是当前节拍陀螺测量角速度。时间积分t(k)＝t(k-1)+ts，t(k)为当前节拍陀螺积分时长，t(k-1)是上一节拍陀螺积分时长，然后跳至步骤6；
步骤5：记录当前星敏感器的姿态，记为θstrf；计算陀螺的常漂残量
计算结束跳至步骤7；
步骤6：判断陀螺积分是否超出阈值，如果是，则跳至步骤8；如果否，则结束；
步骤7：判断陀螺常漂残量Bremain是否小于阈值BC(一般BC不大于星敏感器三轴姿态系统误差的一半，由地面注数决定)，如果是，则跳至步骤9；如果否，则跳至步骤8；
步骤8：清除陀螺修星敏标志flag＝0；(说明：flag是表示是否可以采用陀螺积分值来辨识星敏感器轨道周期系统误差的标志，在第三部分中使用。flag＝0表示陀螺一个轨道周期积分角度太大，不能用于修正星敏感器轨道周期系统误差，反之，flag＝1表示陀螺一个轨道周期积分角度满足要求，可用于修正星敏感器轨道周期系统误差。)计算结束跳至步骤10；
步骤9：建立陀螺修星敏标志flag＝1，计算结束跳至步骤10；
步骤10：计算陀螺常值漂移这里为本次(本轨道周期)计算得到的陀螺常值漂移，其中，ktl∈[0,1]为陀螺常值漂移更新权系数，由地面上注，为上次计算得到的陀螺常值漂移；
积分时间清零t(k)＝0；
陀螺积分清零θgyro(k)＝0；
记录当前轨道钟T0＝torbit；
记录陀螺积分残差θ0＝θgyro(k)-θstrf+θstr0；(说明：记录的陀螺积分残差θ0和记录的当前轨道钟T0＝torbit，均用于第三部分，虽然陀螺一个轨道周期积分角度足够小，满足采用陀螺积分修正星敏感器轨道周期系统误差的条件，但要实现星敏感器系统误差高精度补偿，还需要将陀螺积分残差的影响减至最小。)
第三部分：基于陀螺的星敏感器轨道周期系统误差辨识及补偿表更新
第二部分中，采用星敏感器将陀螺常值漂移修正，并给出陀螺积分残差水平，记录了陀螺积分初始时刻(也是结束时刻，因为陀螺刚好积分一个轨道周期)，给出了是否可以采用陀螺积分辨识星敏感器系统误差的标志。在第三部分中，将给出如何基于上述信息，实现补偿表参数的建立和更新。
同样的道理，与第二部分一致，此处以其中一轴姿态θ为例进行设计说明，如图3所示：
步骤1：判断陀螺修星敏感器的标志flag是否为可用，如果是(flag＝＝1)，则跳至步骤2；如果否，则退出；
步骤2：判断当前轨道钟torbit是否大于记录时刻轨道钟T0，如果是(torbit-T0≥0)，则跳至步骤3；如果否，则跳至步骤4；
步骤3：修正量偏差计算(式中，Δθ为当前节拍基于陀螺积分修正星敏感器时，应扣除的修正量偏差，torbit为当前轨道钟，T0为记录的轨道钟，θ0为陀螺积分残差，T为航天器当前轨道周期。该公式的原理是：由于陀螺一个轨道周期内积分必然存在残差θ0，为了减小以陀螺积分为基准对星敏感器轨道周期系统误差修正的影响，应将此残差θ0的影响消除，其中，(torbit-T0)表示当前点离上个周期陀螺积分结束点的距离，始终为正数，当前轨道钟离陀螺积分结束时刻(也就是记录的轨道钟T0)越近，残差等比例越小，反之，亦然)，计算结束跳至步骤5；
步骤4：修正量偏差计算(也即如果当前轨道钟小于记录的轨道钟，需要进行适应性处理，其中，(T+torbit-T0)表示当前点离上个周期陀螺积分结束点的距离，始终为正数)，计算结束跳至步骤5；
步骤5：计算即将要更新的补偿表行号
当轨道钟与补偿表中轨道钟相等时，就基于最新的星敏感器系统误差信息更新到表格相关位置，在一个完整的轨道周期中，整个补偿表将被更新一轮，每个点仅更新一次。补偿表的行号对应内容更新前，需要计算得到即将要更新的补偿表行号，以便明确本次待更新的是补偿表中的第几行数据。
步骤6：判断当前节拍是否需要对相应的补偿表进行更新，如果满足torbit＝＝(M-1)·ΔT的条件，说明此节拍需要对补偿表进行更新，转至步骤7，否则退出；满足这个条件，说明当前轨道钟时刻与补偿表中的轨道钟(可以理解为补偿表的地址)相对应，此时将识别出的星敏感器系统误差信息更新进入补偿表的对应位置。为后续用插值算法复现星敏常值漂移做准备。经过一整个轨道周期(最后一行与第一行同时更新)，其余补偿表内容均会被更新一轮。
步骤7：计算第M个点的平均修正量：
θxz＝θgyro(k)-Δθ-θstr(k)
其中，θxz是当前节拍计算得到的星敏感器系统误差值，θgyro(k)是当前节拍陀螺积分值，θstr(k)为系统误差经补偿后的星敏感器实测值。
执行第M个点修正参数调整(对第M行数据进行更新)：
其中kxm∈[0,1]，表示修正量更新权系数，由地面上注，横杠表示取平均值，计算结束跳至步骤8；
步骤8：判断当前修正点是否为第一个点，具体判据为：M＝＝1，如果是，则跳至步骤9；如果否，则退出；
步骤9：将第一个点修正值同时存入最后一个点计算结束退出；
轨道钟是周期性的，当轨道钟大小等于轨道周期，轨道钟会清零，而补偿表是线性的，为了将补偿表覆盖整个轨道，补偿表最后一行对应的轨道钟必须大于轨道周期，如果采用和其他点一样的方式，这个点是不可能被更新的。虽然补偿表第一行和最后一行数值上差别很大，但物理上这两行对应的位置是最接近的，因此，在更新补偿表第一行修正值的同时，还需要同时更新补偿表最后一行。
第四部分：基于插值法的星敏感器系统误差补偿算法
第三部分中，得到了星敏感器轨道周期系统误差的补偿表，第四部分将介绍基于补偿表，采用插值算法实现星敏感器轨道周期系统误差的补偿。
同样的道理，与第二部分一致，此处以其中一轴姿态θ为例进行设计说明，如图4所示：
步骤1：根据当前轨道钟torbit，当前轨道钟前一行补偿表内容Block[M]_fy，当前轨道钟后一行补偿表内容Block[M+1]，采用插值算法，计算得到当前轨道钟对应的星敏系统误差：
其中，M的计算方法详见第三部分步骤5。
步骤2：计算得到扣除系统误差后的星敏感器姿态：
θstr(k)＝θstr_o(k)-θbc
其中，θstr_o(k)是未经系统误差补偿的星敏测量值。
以上所有的变量，数据类型不全一致，但其初始值均为0。
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。