[发明专利]一种基于星敏感器的航天器姿态角速度测量方法

摘要

本发明公开了一种基于星敏感器的航天器姿态角速度测量方法，目的是解决现有方法精度受星敏感器测量误差影响较大的问题。技术方案是读入初始时刻t0的星图，获得t0的星敏感器测量矢量及导航星信息；令t＝t0+δt，读入t时刻星图，经星点提取、序列星图识别，对同时出现在前一帧星图和当前帧星图的导航星进行编号，获得同时出现在t-δt和t的导航星总数n以及对应的星敏感器测量矢量对的集合Ω1；然后对卡尔曼滤波器进行初始化，估计航天器姿态角速度，重置滤波器初始值，令t＝t+δt，重复读入t时刻星图，经星点提取、序列星图识别，对导航星进行编号，获得n以及Ω1的步骤。本发明能消除星敏感器测量随机噪声对姿态角速度估计的影响，提高测量的精度。

主权项

1.一种基于星敏感器的航天器姿态角速度测量方法，其特征在于包括以下步骤：第一步，读入初始时刻t₀的星图，进行星点提取后，利用全天球识别模式的星图识别算法获得t₀时刻的星敏感器测量矢量及相应的导航星信息；第二步，令t＝t₀+δt，δt为星敏感器采样时间间隔，读入t时刻星图，进行星点提取、序列星图识别，对同时出现在前一帧星图和当前帧星图的导航星进行编号，获得同时出现在t-δt和t时刻的星图中的导航星总数目n，以及对应的星敏感器测量矢量对的集合Ω1，具体方法如下：2.1读入t时刻星图，进行星点提取，获得t时刻星图中星点的质心坐标，作为星点在星敏感器像面坐标系中的投影位置；2.2进行序列星图识别，获得t时刻星图中识别的导航星i的星敏感器测量矢量v(t)_i，v(t)i=1x2(t)i+y2(t)i+f2x(t)iy(t)if---(1)]]>其中，x(t)_i和y(t)_i分别表示t时刻星图中识别的导航星i在星敏感器像面坐标系中的投影位置的x坐标和y坐标；f表示星敏感器镜头焦距；i为正整数，表示识别的导航星编号；2.3对同时出现在前一帧星图和当前帧星图的导航星进行编号，具体方法为，统计同时出现在前一帧t-δt时刻和当前帧t时刻的星图中的导航星，将总数目记为n，并按赤经由小到大编号为1，2，…k，…n；将导航星k在t-δt星图中的星敏感器测量矢量v(t-δt)_k和t时刻的星图中的星敏感器测量矢量v(t)_k，作为数据对记入集合Ω1，其中k为导航星序号，1≤k≤n；第三步，对卡尔曼滤波器进行初始化：3.1设置估计噪声协方差矩阵初值P₀为3×3的单位矩阵I_3×3；3.2取Ω1中的数据，利用公式(2)，获得t₀+δt时刻的航天器姿态角速度初值ω(t₀+δt)_x、ω(t₀+δt)_y、ω(t₀+δt)_z分别表示t₀+δt时刻的航天器在x、y、x方向上的姿态角速度分量；令X0=ω(t0+δt)xω(t0+δt)yω(t0+δt)z,]]>将X₀作为卡尔曼滤波器的状态矢量的初值；ω(t0+δt)xω(t0+δt)yω(t0+δt)z=-1δt{Σk=1n[v(t0+δk)k×]T[v(t0+δt)k×]}-1Σk=1n[v(t0+δt)k×]Tv(t0)k---(2)]]>式2中，[]^T表示矩阵转置，{}^-1表示矩阵求逆，[v(t₀+δt)_k×]表示t₀+δt时刻导航星k的星敏感器测量矢量v(t0+δt)k=1x2(t0+δt)k+y2(t0+δt)k+f2x(t0+δt)ky(t0+δt)kf]]>的叉乘矩阵，具体形式为：[v(t0+δt)k×]=1x2(t0+δt)k+y2(t0+δt)k+f20-fy(t0+δt)kf0-x(t0+δt)k-y(t0+δt)kx(t0+δt)k0---(3)]]>第四步，估计航天器姿态角速度：4.1取集合Ω1中的数据，以导航星序号k的增加作为时间序列，构建卡尔曼滤波器的状态方程和测量方程如下：X(t)k=I3×3X(t)k-1+wZ(t)k=H(t)kX(t)k+u---(4)]]>其中，X(t)=ω(t)xω(t)yω(t)z]]>表示t时刻的卡尔曼滤波器的状态矢量，ω(t)_x、ω(t)_y、ω(t)_z分别表示t时刻的航天器在x、y、x方向上的姿态角速度分量；X(t)_k和X(t)_k-1分别表示t时刻导航星k和导航星k-1对应的卡尔曼滤波器的状态矢量；I_3×3为3×3单位矩阵；w_k-1表示导航星k-1对应的过程噪声，u_k表示导航星k对应的测量噪声；表示t时刻星敏感器测量矢量k的一阶导；观测矩阵H(t)_k＝[v(t)_k×]，[v(t)_k×]表示t时刻导航星k的星敏感器测量矢量v(t)k=1x2(t)k+y2(t)k+f2x(t)xy(t)yf]]>的叉乘矩阵，具体形式为：[v(t)k×]=1x2(t)k+y2(t)k+f20-fy(t)kf0-x(t)k-y(t)kx(t)k0---(5)]]>4.2输入状态矢量的初值X₀、估计噪声协方差矩阵初值P₀、过程噪声协方差矩阵Q、测量噪声的协方差矩阵R，估计当前时刻的航天器姿态角速度估计值X(t)和估计噪声协方差矩阵P_n，具体步骤如下：4.2.1令k＝1，输入状态矢量的初值X₀、估计噪声协方差矩阵初值P₀；4.2.2采用公式(6)预估导航星k对应的卡尔曼滤波器的状态矢量，记为X(t)k-=I3×3X(t)k-1---(6);]]>4.2.3采用公式(7)预估导航星k对应的噪声协方差矩阵，记为Pk-=Pk-1+Q---(7);]]>4.2.4采用公式(8)计算导航星k对应的卡尔曼增益矩阵，记为K_kKk=Pk-[Pk-+R]-1---(8);]]>4.2.5采用公式(9)更新估计，获得导航星k对应的的卡尔曼滤波器的状态矢量X(t)_kX(t)k=X(t)k-+Kk[Z(t)k-H(t)kX(t)k-]---(9);]]>4.2.6采用公式(10)计算更新后的导航星k对应的估计噪声协方差矩阵P_kPk=[I3×3-KkH(t)k]Pk----(10);]]>4.2.7如果k＝n，清空集合Ω1，输出X(t)_n和P_n，作为当前时刻的航天器姿态角速度估计值X(t)和估计噪声协方差矩阵，转第五步；否则，令k＝k+1，转步骤4.2.2；第五步，重置滤波器初始值：5.1设置估计噪声协方差矩阵初值P₀＝P_n，其中P_n为第四步最后得到的估计噪声协方差矩阵；5.2利用线性插值的方法求下一帧星图对应的航天器姿态角速度预测值，将航天器姿态角速度预测值作为下一帧星图的卡尔曼滤波器的状态矢量的初值X₀，计算公式为：X₀＝2X(t)-X₀    (11)第六步，令t＝t+δt，读入t时刻星图，进行2.1所述的星点提取、2.2所述的序列星图识别，2.3所述的对同时出现在前一帧星图和当前帧星图的导航星进行编号，获得同时出现在t-δt和t时刻的星图中的导航星总数目n，以及对应的星敏感器测量矢量对的集合Ω1，转第四步。