本发明属于遥感卫星地面预处理领域，特别是涉及到一种星敏感器原始观测数据预处理方法及系统。
背景技术：
目前高分辨率光学卫星广泛采用由星敏感器和陀螺构成的卫星姿态确定系统实现精密姿态确定，并用于光学影像的高精度几何处理。星敏感器作为一种常用的高精度姿态测量部件，其主要工作原理是通过恒星相机对恒星进行摄影，利用获取的恒星影像确定卫星的姿态，处理流程包括星点提取、导航星表建立与优化以及星图识别等关键步骤。由于卫星在轨运行过程中所处的空间环境比较复杂以及成像器件本身特性，导致以下干扰因素的产生：(1)星敏感器获取的星图存在各种噪声，例如光子散粒噪声、读出噪声以及暗电流噪声；(2)星图数据纹理信息匮乏、信噪比比较低，存在拖尾现象；(3)星图中导航星分布不均匀；(4)星图匹配识别错误。基于以上因素会导致星敏感器输出的观测值中存在一些不正常观测值，因此需要对星敏感器输出的观测值进行可用性判识，以免将非正常工作中的观测数据引入组合定姿系统中，影响最终精密姿态的确定与光学影像的几何处理。
技术实现要素：
本发明针对星敏感器输出观测值中存在粗差难以探测与修复的问题，提供了一种基于光轴夹角变化检测技术的星敏感器观测数据预处理技术方案。
本发明提供的技术方案为一种星敏感器原始观测数据预处理方法，包括以下步骤：
步骤1，基于陀螺时标信息对多星敏感器原始观测数据的时间基准进行统一；
步骤2，基于步骤1的结果，计算多星敏感器间光轴夹角时间变化序列，多星敏感器间光轴夹角计算实现方式如下，
设在相同时刻t的星敏感器A的四元数观测值为星敏感器B的四元数观测值为得到时刻t的本体到惯性系的星敏感器A、B分别相应旋转矩阵表达式如下，
进一步得到星敏感器A、星敏感器B的光轴矢量ZA,ZB在惯性系中表示如下，
ZA＝[2(q1Aq3A+q2Aq0A) 2(q2Aq3A-q1Aq0A) -q1A2-q2A2+q3A2+q0A2]T
ZB＝[2(q1Bq3B+q2Bq0B) 2(q2Bq3B-q1Bq0B) -q1B2-q2B2+q3B2+q0B2]T
计算t时刻两光轴的夹角αt，
αt＝arccos(ZA·ZB)
步骤3，基于星敏感器地面标定的安装矩阵参数计算星敏感器理论光轴夹角大小，进一步构建星敏感器原始观测数据质量控制判定模型，探测星敏感器中存在的粗差，最后将每个星敏感器的时间序列观测值分成正常输出观测值组与非正常输出观测值组；
所述星敏感器理论光轴夹角计算实现方式如下，
设星敏感器A与星敏感器B的安装矩阵分别表示如下，
则星敏感器A与星敏感器B的光轴在卫星本体坐标系中矢量表示如下，
得到星敏感器A与B的光轴在本体系中的夹角αAB如下，
所述星敏感器原始观测数据的质量控制判定模型如下，
其中，m表示阈值系数，αi表示观测时刻ti的星敏感器A与B的光轴夹角大小，δα为星敏感器A与星敏感器B光轴夹角中误差；
步骤4，根据步骤3的星敏感器所得观测值分组结果，将正常输出观测值作为拟合基准点，基于滑动窗口多项式拟合对非正常输出的星敏感器观测值进行修复。
而且，步骤1中，对多星敏感器原始观测数据的时间基准进行统一采用球面线性插值模型。
而且，步骤3中，误差δα计算如下，
其中，N表示星敏感器观测历元个数。
而且，步骤4中，基于滑动窗口多项式拟合对非正常输出的星敏感器观测值进行修复，实现方式如下，
设某星敏感器A输出的正常输出观测值组包括n个时间序列观测值(q1,q2,q3…qn-1,qn)，tk为非正常输出观测值组中的起始时刻，
记n个历元的姿态四元数为(q0i,q1i,q2i,q3i)i＝1,2,…,n，相应m-1次最佳正交多项式拟合如下，
其中，t表示时刻，aj,j＝0,1,…,m-1表示多项式系数，设上式为各正交多项式δj(t)的线性组合，
其中，cj,j＝0,1,…,m-1表示正交多项式系数；
根据最小二乘原理得tk时刻姿态四元数拟合值如下，
其中，表示tk时刻四元数矢量部分拟合值，表示tk时刻四元数标量部分拟合值，表示四元数矢量部分正交多项式拟合系数，表示四元数矢量部分正交多项式。
本发明相应提供一种星敏感器原始观测数据预处理系统，包括以下模块：
第一模块，用于基于陀螺时标信息对多星敏感器原始观测数据的时间基准进行统一；
第二模块，用于计算多星敏感器间光轴夹角时间变化序列，多星敏感器间光轴夹角计算实现方式如下，
设在相同时刻t的星敏感器A的四元数观测值为星敏感器B的四元数观测值为得到时刻t的本体到惯性系的星敏感器A、B分别相应旋转矩阵表达式如下，
进一步得到星敏感器A、星敏感器B的光轴矢量ZA,ZB在惯性系中表示如下，
ZA＝[2(q1Aq3A+q2Aq0A) 2(q2Aq3A-q1Aq0A) -q1A2-q2A2+q3A2+q0A2]T
ZB＝[2(q1Bq3B+q2Bq0B) 2(q2Bq3B-q1Bq0B) -q1B2-q2B2+q3B2+q0B2]T
计算t时刻两光轴的夹角αt，
αt＝arccos(ZA·ZB)
第三模块，用于基于星敏感器地面标定的安装矩阵参数计算星敏感器理论光轴夹角大小，进一步构建星敏感器原始观测数据质量控制判定模型，探测星敏感器中存在的粗差，最后将每个星敏感器的时间序列观测值分成正常输出观测值组与非正常输出观测值组；
所述星敏感器理论光轴夹角计算实现方式如下，
设星敏感器A与星敏感器B的安装矩阵分别表示如下，
则星敏感器A与星敏感器B的光轴在卫星本体坐标系中矢量表示如下，
得到星敏感器A与B的光轴在本体系中的夹角αAB如下，
所述星敏感器原始观测数据的质量控制判定模型如下，
其中，m表示阈值系数，αi表示观测时刻ti的星敏感器A与B的光轴夹角大小，δα为星敏感器A与星敏感器B光轴夹角中误差；
第四模块，用于根据第三模块的星敏感器所得观测值分组结果，将正常输出观测值作为拟合基准点，基于滑动窗口多项式拟合对非正常输出的星敏感器观测值进行修复。
而且，第一模块中，对多星敏感器原始观测数据的时间基准进行统一采用球面线性插值模型。
而且，第三模块中，误差δα计算如下，
其中，N表示星敏感器观测历元个数。
而且，第四模块中，基于滑动窗口多项式拟合对非正常输出的星敏感器观测值进行修复，实现方式如下，
设某星敏感器A输出的正常输出观测值组包括n个时间序列观测值(q1,q2,q3…qn-1,qn)，tk为非正常输出观测值组中的起始时刻，
记n个历元的姿态四元数为(q0i,q1i,q2i,q3i)i＝1,2,…,n，相应m-1次最佳正交多项式拟合如下，
其中，t表示时刻，aj,j＝0,1,…,m-1表示多项式系数，设上式为各正交多项式δj(t)的线性组合，
其中，cj,j＝0,1,…,m-1表示正交多项式系数；
根据最小二乘原理得tk时刻姿态四元数拟合值如下，
其中，表示tk时刻四元数矢量部分拟合值，表示tk时刻四元数标量部分拟合值，表示四元数矢量部分正交多项式拟合系数，表示四元数矢量部分正交多项式。
本发明提供了一种基于光轴夹角变化检测技术的星敏感器原始观测数据预处理技术方案，实现了星敏感器非正常输出观测值的有效辨识与修复，通过本发明提供的技术方案可以有效削弱星敏感器粗差观测值影响，进一步实现高精度姿态确定，为高分辨率光学影像高精度几何处理提供保证。
附图说明
图1为本发明实施例的流程图。
具体实施方式
以下结合附图和实施例详细说明本发明技术方案。
图1表示的是基于光轴夹角变化检测技术的星敏感器原始观测数据预处理方法实施流程图，以下针对实施例流程中的各步骤，对本发明方法做进一步详细描述。
步骤1，基于陀螺输出观测信息，对多星敏感器观测值的时间基准进行统一，保证多星敏感器输出观测值时刻一致性。
通常包括两个或以上星敏感器，处理方式类似。陀螺输出观测信息即陀螺输出观测数据，多星敏感器观测值即多星敏感器输出观测值。基于这两者实现对多星敏感器观测值的时间基准进行统一，得到星敏感器时间同步数据。陀螺作为星体的短期姿态参考，能够连续提供星体的三轴姿态角速度信息，但由于多台星敏感器与陀螺并不能实现完全时间同步，为了后续在步骤2得到多星敏感器间光轴夹角变化时间序列，本发明基于陀螺的时标信息采用球面线性插值模型实现多星敏感器间原始观测数据时间同步。具体实现方式如下：
球面线性插值不仅会保持其单位长度不变，并且若将q1和q2看为四维空间单位球上的两个点，其将以恒定的速率扫过q1和q2之间的夹角θ。假设t时刻所需插值计算的四元数q(t)位于连接q1和q2的弧上,与q1构成的夹角为tθ(t∈[0,1])与q2构成的夹角为(1-t)θ(t∈[0,1]),于是q(t)可表示为
q(t)＝C1(t)q1+C2(t)q2 (1)
根据相似三角形关系,可以推导出系数C1(t)和C2(t)的表达式：
即：
步骤2，计算多星敏感器间光轴夹角时间序列数据：基于步骤1的结果，进一步计算多星敏感器间光轴夹角时间变化序列，得到星敏感器间光轴夹角数据。当星敏观测误差满足正态分布，那么时间序列光轴夹角变化也满足正态分布，夹角期望为无偏估计。多星敏感器间光轴夹角数据计算方法的具体实现方式如下：
假设在相同时刻t的星敏感器A的四元数观测值为星敏感器B的四元数观测值为得到时刻t的本体到惯性系的星敏感器A、B分别相应旋转矩阵表达式如下：
进一步得到星敏感器A、星敏感器B的光轴矢量ZA,ZB在惯性系中表示如下：
ZA＝[2(q1Aq3A+q2Aq0A) 2(q2Aq3A-q1Aq0A) -q1A2-q2A2+q3A2+q0A2]T (5)
ZB＝[2(q1Bq3B+q2Bq0B) 2(q2Bq3B-q1Bq0B) -q1B2-q2B2+q3B2+q0B2]T
可计算t时刻两光轴的夹角αt：
αt＝arccos(ZA·ZB) (6)
故基于上述算法，对于时间变化中每个时刻分别计算光轴间的夹角，可以得到多星敏感器间光轴夹角变化时间序列，为步骤3实施的前提条件。
步骤3，星敏感器观测值粗差探测：基于星敏感器地面标定的安装矩阵参数解算星敏感器光轴夹角大小，进一步构建星敏感器原始观测数据质量控制判定模型，探测星敏感器中可能存在的粗差，最终将每个星敏感器的时间序列观测值分成正常输出观测值组与非正常输出观测值组，得到分组后的星敏感器观测值序列。该步骤包括以下子步骤：
步骤3.1，基于地面标定的星敏感器的安装矩阵参数计算星敏感器之间的光轴夹角大小，具体实现方式如下：
设星敏感器A与星敏感器B的安装矩阵分别表示如下：
则星敏感器A与星敏感器B的光轴在卫星本体坐标系中矢量可以表示成：
同样可以得到星敏感器A与B的光轴在本体系中的夹角αAB：
步骤3.2，计算星敏感器A与星敏感器B光轴夹角中误差δα，具体实现方式如下：
其中，αi表示观测时刻ti的星敏感器A与B的光轴夹角大小，通过步骤2得到；N表示星敏感器观测历元个数。
步骤3.3，构建星敏感器观测数据质量控制判定模型，探测星敏感器中可能存在的粗差，最终将每个星敏感器的时间序列观测值分成正常输出观测值组与非正常输出观测值组。具体实现方式如下：
星敏感器原始观测数据质量控制判定模型如下：
其中，m表示阈值系数，且1≤m≤3。故依据上述构建的判定模型对星敏感器输出的观测值时间序列进行判定，将正常与非正常输出的观测值各自构成一组时间序列，即分组后的星敏感器观测值序列。
步骤4，星敏感器观测值数据修复，得到多星敏感器可靠观测值数据：根据步骤3的星敏感器A与B时间序列观测值分组结果，将正常输出观测值作为拟合基准点，基于滑动窗口多项式拟合对非正常输出的星敏感器原始观测值进行修复。具体过程如下：
以星敏感器A为例，设其输出的正常观测值序列组包括n个时间序列观测值(q1,q2,q3…qn-1,qn)，非正常输出观测值的时刻序列组根据实际分组情况得到，例如(tk,tk+2,tk+6…,tk+l)，tk为非正常输出观测值时刻序列组中的起始时刻，采用的拟合多项式模型如下：
记n个历元的姿态四元数为(q0i,q1i,q2i,q3i)i＝1,2,…,n，其中qri,r＝1,2,3参数的m-1次最佳正交多项式拟合为：
其中，t表示时刻，aj,j＝0,1,…,m-1表示多项式系数。设上式为各正交多项式δj(t)(j＝0,1,…,m-1)的线性组合：
其中，cj,j＝0,1,…,m-1表示正交多项式系数。
其中δj(t)可用递推公式构造：
δ0(t)＝1；δ1(t)＝(t-α1)；
δj(t)＝(t-αj)δj-1(t)-βjδj-2(t) (14)
j＝2,…,m-1
根据最小二乘原理可得tk时刻姿态四元数拟合值：
其中，dj表示正交多项式平方和，表示tk时刻四元数矢量部分拟合值，表示tk时刻四元数标量部分拟合值，表示四元数矢量部分正交多项式拟合系数，表示四元数矢量部分正交多项式。故根据上述模型对非正常输出观测值的时刻序列进行修复。
具体实施时，本发明所提供方法可基于软件技术实现自动运行流程，也可采用模块化方式实现相应系统。
本发明相应提供一种星敏感器原始观测数据预处理系统，包括以下模块：
第一模块，用于基于陀螺时标信息对多星敏感器原始观测数据的时间基准进行统一；
第二模块，用于计算多星敏感器间光轴夹角时间变化序列，多星敏感器间光轴夹角计算实现方式如下，
设在相同时刻t的星敏感器A的四元数观测值为星敏感器B的四元数观测值为得到时刻t的本体到惯性系的星敏感器A、B分别相应旋转矩阵表达式如下，
进一步得到星敏感器A、星敏感器B的光轴矢量ZA,ZB在惯性系中表示如下，
ZA＝[2(q1Aq3A+q2Aq0A) 2(q2Aq3A-q1Aq0A) -q1A2-q2A2+q3A2+q0A2]T
ZB＝[2(q1Bq3B+q2Bq0B) 2(q2Bq3B-q1Bq0B) -q1B2-q2B2+q3B2+q0B2]T
计算t时刻两光轴的夹角αt，
αt＝arccos(ZA·ZB)
第三模块，用于基于星敏感器地面标定的安装矩阵参数计算星敏感器理论光轴夹角大小，进一步构建星敏感器原始观测数据质量控制判定模型，探测星敏感器中存在的粗差，最后将每个星敏感器的时间序列观测值分成正常输出观测值组与非正常输出观测值组；
所述星敏感器理论光轴夹角计算实现方式如下，
设星敏感器A与星敏感器B的安装矩阵分别表示如下，
则星敏感器A与星敏感器B的光轴在卫星本体坐标系中矢量表示如下，
得到星敏感器A与B的光轴在本体系中的夹角αAB如下，
所述星敏感器原始观测数据的质量控制判定模型如下，
其中，m表示阈值系数，αi表示观测时刻ti的星敏感器A与B的光轴夹角大小，δα为星敏感器A与星敏感器B光轴夹角中误差；
第四模块，用于根据第三模块的星敏感器所得观测值分组结果，将正常输出观测值作为拟合基准点，基于滑动窗口多项式拟合对非正常输出的星敏感器观测值进行修复。
各模块具体实现可参见相应步骤，本发明不予赘述。
本文中所描述的具体实例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。