CN107449444B - 一种多星图姿态关联的星敏感器内参数标定方法 - Google Patents

一种多星图姿态关联的星敏感器内参数标定方法 Download PDF

Info

Publication number
CN107449444B
CN107449444B CN201710581215.0A CN201710581215A CN107449444B CN 107449444 B CN107449444 B CN 107449444B CN 201710581215 A CN201710581215 A CN 201710581215A CN 107449444 B CN107449444 B CN 107449444B
Authority
CN
China
Prior art keywords
star
star sensor
sensor
calibration
coordinate system
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710581215.0A
Other languages
English (en)
Other versions
CN107449444A (zh
Inventor
戴东凯
刘宇
马丽衡
秦石乔
王省书
吴伟
郑佳兴
黄宗升
魏文俭
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN201710581215.0A priority Critical patent/CN107449444B/zh
Publication of CN107449444A publication Critical patent/CN107449444A/zh
Application granted granted Critical
Publication of CN107449444B publication Critical patent/CN107449444B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C25/00—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass

Landscapes

  • Engineering & Computer Science (AREA)
  • Manufacturing & Machinery (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Navigation (AREA)

Abstract

一种多星图姿态关联的星敏感器内参数标定方法及其装置，将GPS天线3、陀螺组合体1、待标定星敏感器2分别与GPS接收机4通信，将数据处理计算机5分别与所述陀螺组合体1、待标定星敏感器2连接；利用GPS接收机4采集得到的世界协调时（UTC）时间信息实现陀螺组合体1和星敏感器2数据同步采集并在数据处理计算机5中完成待标定星敏感器2的标定算法求解。该方法实现了利用陀螺组合体提供精确旋转角度信息实现多帧星图的拼接关联，从而增加用于星敏感器内参数标定的观测数据样本，可以用于动态条件下的标定，并对待标定星敏感器2的运动状态无严格的要求，提高了标定精度和可靠性，简化了标定流程且便于实施待标定星敏感器2的动态标定，该方法和装置也适用于惯性/天文组合导航系统的标定。

Description

一种多星图姿态关联的星敏感器内参数标定方法
技术领域
本发明涉及航天测量领域中的星敏感器内参数标定方法及其装置；具体地说是一种姿态关联帧星敏感器的内参数标定方法及其专用装置。
背景技术
星敏感器是航天器姿态控制的核心传感器之一，它通过敏感星光的入射角度信息实现航天器姿态的测量。对光学镜头焦距、主点以及畸变等参数(统称为内参数)的精确标定是星敏感器实现高精度姿态测量的前提和关键，因此星敏感器在投入使用前都必须在地面进行严格的标定。
我国的国家军用标准“星敏感器标定与测试方法”(GJB 8237-2013)规定了星敏感器内参数标定的标准方法，其实施过程需要在光学实验室条件下采用单星模拟器配合高精度二维转台实现标定，标定精度高，但是对标定环境和设备的要求苛刻，成本较高。另一种通用的标定方法是根据恒星对间的星角距正交变换不变的特性，以采集到的星点的图像坐标和对应观测时刻导航星的视赤经和视赤纬作为输入，无需高精度标定设备即可实现星敏感器镜头参数的估计。然而，这种方法没有充分考虑星角距误差的真实统计特性，简单地把它当作白噪声进行处理，因而算法在理论上不是最优估计，标定精度较差。
如何提高星敏感器内参数的标定精度、简化标定流程从而实现星敏感器姿态测量精度整体提升是本领域技术人员极为关注的技术问题。
发明内容
本发明在国军标GJB 8237-2013规定的星敏感器标定数学模型基础上进行改进，在大气扰动小的地点对晴朗夜空进行恒星拍摄，利用陀螺组合体(Gyroscope Unit,GU)提供精确旋转角度信息实现多帧星图的拼接关联，从而增加用于星敏感器内参数标定的观测数据样本，可以用于动态条件下的标定，提高标定精度和可靠性，简化了标定流程且便于实施。
为解决上述问题，本发明所采用的技术方案是：一种利用多星图姿态关联实现星敏感器内参数标定方法采用的装置由陀螺组合体、GPS天线、GPS接收机、数据处理计算机以及待标定星敏感器组成；
所述陀螺组合体与所述待标定星敏感器刚性安装，以下将刚性安装的陀螺组合体和星敏感器简称为星惯组合系统；
所述GPS天线与所述GPS接收机通信；
所述陀螺组合体与所述GPS接收机通信；
所述待测星敏感器与所述GPS接收机通信，利用GPS接收机采集得到的世界协调时(UTC)时间信息可以实现陀螺组合体和星敏感器数据的时间同步；
所述数据处理计算机分别与所述陀螺组合体、待测星敏感器连接，并进行同步数据采集，在数据处理计算机中完成星敏感器的标定算法求解。
为实现多星图姿态关联的星敏感器内参数标定方法的技术方案采用以下步骤实现：
标定前选择在晴朗无云的夜晚、远离城市灯光的环境下实施标定，将星惯组合系统平放于地面，使星敏感器的光轴大致朝向天顶方向，确保星敏感器的视场内无障碍物遮挡以便于拍摄恒星，标定开始前将陀螺组合体通电1小时以上；
步骤1：标定开始；
1.1通过数据处理计算机控制陀螺组合体，使其开始进行姿态解算，并输出陀螺组合体相对于惯性坐标系的姿态，利用数据处理计算机开始采集陀螺组合体输出的姿态数据和每个数据所对应的UTC时间信息，陀螺组合体进行姿态解算的算法如下：
1.1.1给定初始时刻t0陀螺组合体的初始姿态矩阵为
Figure GDA0002340045830000021
1.1.2第k个采样时刻记为tk，陀螺组合体中三个陀螺11，陀螺12，陀螺13输出的角增量分别记为
Figure GDA0002340045830000022
利用陀螺的角增量构造姿态四元数：
Figure GDA0002340045830000023
其中
Figure GDA0002340045830000024
符号[]T表示矩阵的转置；令q0、q1、q2、q3分别为qk的第1～4个元素；利用式(1)初始姿态矩阵计算tk-1时刻到tk时刻的陀螺组合体坐标系的姿态变化矩阵：
Figure GDA0002340045830000025
利用步骤1.1.1给出的初始姿态和式(2)给出的姿态变化矩阵
Figure GDA0002340045830000026
进行陀螺组合体姿态迭代更新，直至测量结束；姿态更新迭代方程如下：
Figure GDA0002340045830000027
(1.2)利用数据处理计算机控制星敏感器使其开始拍摄恒星，并采集星敏感器拍摄得到的星图，记录每帧星图对应的UTC时间。
步骤2：标定数据采集；
在标定过程中通过手动或者二轴转台缓慢地绕任意3个不同方向的转轴分别转动星惯组合系统，使星敏感器在不同的姿态下拍摄恒星，转动星惯组合系统的过程保证星敏感器的光轴的俯仰角不小于45°，在该过程中采用数据处理计算机连续地同步采集星敏感器和陀螺组合体的测量数据，标定数据采集的过程持续小于5分钟。
步骤3：结束标定数据采集；
通过数据处理计算机控制陀螺组合体，使其停止测量，控制星敏感器使其停止拍摄恒星。
步骤4：标定数据处理；
通过数据处理计算机运行以下步骤算法进行标定数据处理，标定数据处理算法如下：
步骤4.1标定数据的预处理；
4.1.1计算载体的平均角速度；步骤2所采集到的陀螺数据相对应的时标记为tk，利用陀螺输出的角增量计算tk时刻的载体运动的近似平均角速度，计算方法如下：
Figure GDA0002340045830000031
其中(m＝x,y,z)；N为1秒内陀螺的采样个数；
Figure GDA0002340045830000032
4.1.2时标取齐；给定步骤2拍摄得到的第i帧星图(i＝1…L，L为全部星图帧数)，相对应的时标记为ti，读取全部陀螺组合体数据时标，查找陀螺组合体时标数据，如果满足|ti-tk|＜τ，(τ为陀螺数据采样间隔)，则保存新的数据记录，该新数据记录的内容包括时标ti，第i帧星图，tk时刻对应的陀螺组合体姿态输出
Figure GDA0002340045830000033
以及载体角速度ωk；
4.1.3剔除大动态数据；
依次读取步骤4.1.2保存的新数据记录，如果ωk＞0.1°/s，则剔除该数据记录；
步骤4.2星点提取与全天时星图识别；
对步骤4.1保存的数据记录中的星图依次进行星点提取和星图识别，星点提取的方法参照《光学技术》2009年第35卷第3期刊载的“基于背景自适应预测的星点提取算法”，提取得到的第j帧的中的第k个星点图像坐标记为(uj,k,vj,k)，其中k＝1…Mj，Mj为第j帧星图中的星点总数；星图识别方法参照《光学精密工程》2009年第17卷第1期刊载的“改进的基于主星的星图识别算法”，通过星图识别可以得到第j帧的中的第k个星点的天球坐标为
Figure GDA0002340045830000034
步骤4.3给定星敏感器星的初始参数矢量
Figure GDA0002340045830000035
和初始畸变参数
Figure GDA0002340045830000036
根据星敏感器光学镜头的标称参数给定星敏感器主点、焦距的初值
Figure GDA0002340045830000037
Figure GDA0002340045830000038
不考虑畸变参数，根据实际星敏感器与陀螺组合体的安装角度关系给定星敏坐标系S相对于陀螺组合体的安装角初值
Figure GDA0002340045830000039
利用给定的星敏感器内参数初步计算第1帧星图对应的星敏坐标系S姿态欧拉角初值
Figure GDA00023400458300000310
令初始参数矢量为
Figure GDA00023400458300000311
给定初始畸变参数
Figure GDA00023400458300000312
步骤4中4.4估计星敏感器的主点、焦距以及外参数；
4.4.1建立姿态帧关联星点成像模型；
建立星敏感器的星敏坐标系S星点成像模型如下：
Figure GDA0002340045830000041
其中，(uj,k,vj,k)为第j帧星图中第k个星点的图像坐标，(u0,v0)、f分别为待标定的星敏感器主点、焦距；
Figure GDA0002340045830000042
为第j帧的中的第k个星点在星敏坐标系S下的坐标；光学镜头畸变采用下式(6)给出的参数模型计算；
Figure GDA0002340045830000043
其中，p1,p2,q1,q2,q3为待标定的星敏感器的镜头畸变参数，令畸变参数矢量为
Kd＝[p1,p2,q1,q2,q3]T，rj,k为星点图像坐标(uj,k,vj,k)相对于主点(u0,v0)的距离；
对于第1帧星图，星点在星敏坐标系S与惯性坐标系下的坐标存在如下映射关系：
Figure GDA0002340045830000044
其中，
Figure GDA0002340045830000045
为第1帧星图中的第k个星点在惯性坐标系i下的坐标；
Figure GDA0002340045830000046
为第1帧星图时刻星敏坐标系S相对于惯性坐标系i的姿态矩阵，可以用欧拉角
Figure GDA0002340045830000047
表示如下：
Figure GDA0002340045830000048
第j帧的中的第k个星点在星敏坐标系S下的坐标可以通过第1帧星图的坐标递推得到：
Figure GDA0002340045830000049
其中，
Figure GDA00023400458300000410
为第j帧星图中的第k个星点在惯性坐标系下的坐标；
Figure GDA00023400458300000411
表示星敏坐标系S相对于陀螺组合体的安装矩阵，可以用欧拉角θ1,θ2,θ3表示如下：
Figure GDA00023400458300000412
Figure GDA00023400458300000413
为陀螺组合体坐标系b从t1时刻到tj时刻的姿态变换矩阵，可以由步骤1.1所述陀螺组合体解算得到t1和tj采样时刻陀螺组合体坐标系b相对于惯性坐标系i的姿态
Figure GDA00023400458300000414
和
Figure GDA00023400458300000415
根据式(9)计算：
Figure GDA00023400458300000416
联合上式可以建立以下成像模型：
Figure GDA0002340045830000051
Figure GDA0002340045830000052
和
Figure GDA0002340045830000053
分别为uj,k、vj,k关于参数u0,v0,f,θ1,θ2,θ3,
Figure GDA0002340045830000054
的函数；
4.4.2由步骤4.3给定星敏感器参数矢量初值
Figure GDA0002340045830000055
畸变参数矢量初值
Figure GDA0002340045830000056
和星点在惯性坐标系i下的坐标
Figure GDA0002340045830000057
根据式上式(10)给出的星点成像模型计算全部星点的图像坐标估计值，根据提取得到的星点图像坐标(uj,k,vj,k)计算星点图像坐标估计误差：
Figure GDA0002340045830000058
Figure GDA0002340045830000059
为利用式(10)计算得到的第j帧星图中第k个星点的图像坐标估计值，
Figure GDA00023400458300000510
为第j帧星图中第k个星点的图像坐标估计误差。
4.4.3估计星敏感器的内外参数误差；
不考虑镜头畸变，星敏感器的参数误差矢量
Figure GDA00023400458300000511
与星点图像坐标估计误差
Figure GDA00023400458300000512
满足如下关系：
Figure GDA00023400458300000513
其中，
Figure GDA00023400458300000514
为星敏感器的参数矢量；
Figure GDA00023400458300000515
为X的估计误差矢量；
Figure GDA00023400458300000516
为
Figure GDA00023400458300000517
对参数矢量X的偏导数矢量；
Figure GDA00023400458300000518
为
Figure GDA00023400458300000519
对参数矢量X的偏导数矢量；
对于步骤(4)中4.2提取的所有的星图的全部星点联立方程，完成步骤(4)中(4.3)；
Figure GDA00023400458300000520
其中：
Figure GDA0002340045830000061
N为全部星图的总帧数，M为第N帧星图的星点数；
将全部星点联立方程简写为Z＝H·ΔX，其中
Figure GDA0002340045830000062
采用最小二乘算法计算ΔX如下：
ΔX＝(HT·H)-1HT·Z (14)
4.4.4更新星敏感器的星敏感器内外参数，完成步骤(4)中(4.4)；
利用估计得到的星敏感器内外参数误差更新星敏感器的内外参数，更新方法如下：
Figure GDA0002340045830000063
然后将新的星敏感器的参数矢量X赋值给内外参数初值
Figure GDA0002340045830000064
即
Figure GDA0002340045830000065
步骤4.5星敏感器的光学镜头畸变估计；
4.5.1给定星敏感器参数矢量
Figure GDA0002340045830000066
畸变参数矢量
Figure GDA0002340045830000067
和星点在惯性坐标系i下的坐标
Figure GDA0002340045830000068
根据式(10)给出的星点成像模型计算全部星点的图像坐标估计值
Figure GDA0002340045830000069
根据步骤4.2提取得到的星点图像坐标(uj,k,vj,k)计算星点图像坐标估计误差，完成步骤(4)中(4.5)：
Figure GDA00023400458300000610
4.5.2估计星敏感器镜头畸变，优化内外参数；
星点图像坐标估计误差
Figure GDA00023400458300000611
满足如下关系：
Figure GDA00023400458300000612
其中，
Figure GDA0002340045830000071
为
Figure GDA0002340045830000072
对参数矢量X的偏导数矢量；
Figure GDA0002340045830000073
为
Figure GDA0002340045830000074
对参数矢量X的偏导数矢量，
Figure GDA0002340045830000075
ΔKd为畸变参数误差；
对于步骤4.2提取到的所有的星图的全部星点联立方程：
Figure GDA0002340045830000076
其中：
Figure GDA0002340045830000077
N为全部星图的总帧数，M为第N帧星图的星点数；
将上式全部星点联立方程简写为Z＝H·ΔX+，其中
Figure GDA0002340045830000078
采用最小二乘算法计算ΔX如下：
ΔX＝(HT·H)-1HT·Z (19)
4.5.3更新星敏感器内外参数和畸变参数；
利用估计得到的星敏感器内外参数误差更新星敏感器的内外参数，更新方法如下：
Figure GDA0002340045830000079
然后将新的星敏感器的参数矢量X赋值给内外参数初值
Figure GDA00023400458300000710
即
Figure GDA00023400458300000711
将新的星敏感器的畸变参数矢量Kd赋值给畸变参数初值
Figure GDA00023400458300000712
步骤4.6判断是否满足标定精度的要求；
利用标定得到的星敏感器参数矢量X、畸变参数矢量Kd和星点在惯性坐标系i下的坐标
Figure GDA00023400458300000713
根据式(10)给出的星点成像模型计算全部星点的图像坐标估计值
Figure GDA0002340045830000081
根据提取得到的星点图像坐标(uj,k,vj,k)计算星点图像坐标估计误差，，完成步骤(4)中(4.6)：
Figure GDA0002340045830000082
如果
Figure GDA0002340045830000083
则跳转到步骤5.7，否则跳转到步骤4.4.2重新进行参数估计；其中max()为取最大值函数，||用于求矢量的模，ε为星点重构误差的阈值。
步骤4.7结束标定。
本发明具有以下技术效果：
1.本发明所述的标定方法无需昂贵的内场标定设备，只需要通过室外拍星便可以完成标定，可以降低标定的硬件成本；
2.该方法可以实现星敏感器的动态标定，对星敏感器的运动状态无严格的要求，可以简化标定步骤，提高标定效率，还可以推广用于星敏感器的在轨标定；
3.可以同时标定星敏感器相对于陀螺组合体的安装角，特别适用于惯性/天文组合导航系统的标定。
附图说明
图1是本发明所述的星敏感器标定装置结构示意图；
图中：1.陀螺组合体，2.星敏感器，3.GPS天线，4.GPS接收机，5.数据处理计算机；
图2是本发明所述陀螺组合体的示意图，陀螺组合体包括三个陀螺：陀螺11，陀螺12，陀螺13；
图3是本发明所述的星敏感器标定步骤；
图4为一组用于星敏感器标定的典型的星惯组合系统转动次序；
图5是本发明所述的星敏感器标定数据处理算法流程图示意图。
具体实施方式
下面将结合附图和实例对本发明作进一步的详细说明；
本发明采用的星敏感器内参数标定装置如图1所示，由陀螺组合体1、待标定星敏感器2、GPS天线3、GPS接收机4、数据处理计算机5组成。所述陀螺组合体1与所述待标定星敏感器2刚性安装，以下将刚性安装的陀螺组合体1和星敏感器2简称为星惯组合系统；所述GPS天线3与所述GPS接收机4通信；所述陀螺组合体1与所述GPS接收机4通信；所述待测星敏感器2与所述GPS接收机4通信，利用GPS接收机4采集得到的世界协调时(UTC)时间信息可以实现陀螺组合体1和星敏感器2数据的时间同步；所述数据处理计算机5分别与所述陀螺组合体1、星敏感器2连接，并进行同步数据采集，在数据处理计算机5中完成星敏感器2的标定算法求解。
如图2所示，所述的陀螺组合体1由三个正交安装的陀螺11、12、13组成，优选的陀螺类型是激光陀螺；
定义陀螺组合体1的坐标系为b系Ob-xbybzb，定义星敏感器2的坐标系为s系Os-xsyszs，惯性坐标系为i系Oi-xiyizi；
图3为本发明测量方法总体流程图，本发明具体实施步骤如下：
标定前选择在晴朗无云的夜晚、远离城市灯光的环境下实施标定；将星惯组合系统平放于地面，使星敏感器2的光轴大致朝向天顶方向，确保星敏感器2的视场内无障碍物遮挡以便于拍摄恒星；标定开始前将陀螺组合体1通电1小时以上。
步骤1：标定开始；
1.1通过数据处理计算机控制控制陀螺组合体1，使其开始进行姿态解算，并输出陀螺组合体1相对于惯性坐标系的姿态，利用数据处理计算机5开始采集陀螺组合体1输出的姿态数据
Figure GDA0002340045830000093
和每个数据所对应的UTC时间信息；
1.2控制星敏感器2使其开始拍摄恒星，利用数据处理计算机5采集待标定星敏感器2拍摄得到的星图，并记录每帧星图对应的UTC时间；
步骤2：标定数据采集；在标定过程中通过手动或者二轴转台缓慢地绕任意3个不同方向的转轴分别转动星惯组合系统，使星敏感器2在不同的姿态下拍摄恒星，转动星惯组合系统的过程保证星敏感器2的光轴的俯仰角不小于45°，在该过程中采用数据处理计算机5连续地同步采集星敏感器2和陀螺组合体1的测量数据；标定数据采集的过程持续小于5分钟；
图4为为一组用于星敏感器标定的典型的星惯组合系统转动次序，可以看到星惯组合系统分别围绕独立的三个正交轴X轴、Y轴以及Z轴进行了三次旋转，且旋转的角度相同；
步骤3：结束标定数据采集；通过数据处理计算机5控制陀螺组合体1，使其停止测量，控制星敏感器1使其停止拍摄恒星；
步骤4：标定数据处理；通过数据处理计算机5运行以下算法进行标定数据处理；标定数据处理算法的流程图如图5所示，具体说明如下：
4.1标定数据的预处理；
4.1.1计算载体的平均角速度；步骤2所采集到的陀螺数据相对应的时标记为tk，利用陀螺输出的角增量计算tk时刻的载体运动的近似平均角速度，计算方法如下：
Figure GDA0002340045830000091
其中(m＝x,y,z)，N为1秒内陀螺的采样个数，
Figure GDA0002340045830000092
4.1.2时标取齐；给定步骤2拍摄得到的第i帧星图(i＝1…L，L为全部星图帧数)，相对应的时标记为ti，读取全部陀螺组合体1数据时标，查找陀螺组合体1时标数据，如果满足|ti-tk|＜τ，(τ为陀螺数据采样间隔)，则保存新的数据记录，该新数据记录的内容包括时标ti，第i帧星图，tk时刻对应的陀螺组合体1姿态输出
Figure GDA0002340045830000101
以及载体角速度ωk；
4.1.3剔除大动态数据；
依次读取步骤4.1.2保存的新数据记录，如果ωk＞0.1°/s，则剔除该数据记录；
4.2星点提取与全天时星图识别；
对步骤4.1保存的数据记录中的星图依次进行星点提取和星图识别，星点提取的方法参照《光学技术》2009年第35卷第3期刊载的“基于背景自适应预测的星点提取算法”，提取得到的第j帧的中的第k个星点图像坐标记为(uj,k,vj,k)，其中k＝1...Mj，Mj为第j帧星图中的星点总数；星图识别方法参照《光学精密工程》2009年第17卷第1期刊载的“改进的基于主星的星图识别算法”，通过星图识别可以得到第j帧的中的第k个星点的天球坐标为
Figure GDA0002340045830000102
4.3给定星敏感器星的主点、焦距以及外参数初值；根据星敏感器2光学镜头的标称参数给定星敏感器2主点、焦距的初值
Figure GDA0002340045830000103
Figure GDA0002340045830000104
不考虑畸变参数，根据实际星敏感器2与陀螺组合体1的安装角度关系给定星敏感器2相对于陀螺组合体1的安装角初值
Figure GDA0002340045830000105
利用给定的星敏感器2内参数初步计算第1帧星图对应的星敏感器2姿态欧拉角初值
Figure GDA0002340045830000106
令初始参数矢量为
Figure GDA0002340045830000107
给定初始畸变参数
Figure GDA0002340045830000108
4.4估计星敏感器2的主点、焦距以及外参数；
4.4.1建立姿态帧关联星点成像模型；
建立星敏感器2的星敏坐标系S星点成像模型如下：
Figure GDA0002340045830000109
其中，(uj,k,vj,k)为第j帧星图中第k个星点的图像坐标，(u0,v0)、f分别为待标定星敏感器2的星敏坐标系S主点、焦距；
Figure GDA00023400458300001010
为第j帧的中的第k个星点在星敏坐标系S下的坐标；光学镜头畸变采用下式(6)给出的参数模型计算；
Figure GDA0002340045830000111
其中p1,p2,q1,q2,q3为待标定的星敏感器2的星敏坐标系S镜头畸变参数，令畸变参数矢量为Kd＝[p1,p2,q1,q2,q3]T，rj,k为星点图像坐标(uj,k,vj,k)相对于主点(u0,v0)的距离；
对于第1帧星图，星点在星敏感器2的星敏坐标系S坐标系与惯性坐标系i下的坐标存在如下映射关系：
Figure GDA0002340045830000112
其中，
Figure GDA0002340045830000113
为第1帧星图中的第k个星点在惯性坐标系i下到坐标；
Figure GDA0002340045830000114
为第一帧星图时刻星敏感器2的星敏坐标系S坐标系相对于惯性坐标系i的姿态矩阵，可以用欧拉角
Figure GDA0002340045830000115
表示如下：
Figure GDA0002340045830000116
第j帧的中的第k个星点在星敏坐标系S下的坐标可以通过第1帧星图的坐标递推得到：
Figure GDA0002340045830000117
其中，
Figure GDA0002340045830000118
为第j帧星图中的第k个星点在惯性坐标系i下到坐标；
Figure GDA0002340045830000119
表示星敏感器2的星敏坐标系S相对于陀螺组合体1的安装矩阵，可以用欧拉角θ1,θ2,θ3表示如下：
Figure GDA00023400458300001110
Figure GDA00023400458300001111
为陀螺组合体坐标系b从t1时刻到tj时刻的姿态变换矩阵，可以由步骤1.1所述陀螺组合体解算得到t1和tj采样时刻陀螺组合体坐标系b相对于惯性坐标系i的姿态
Figure GDA00023400458300001112
和
Figure GDA00023400458300001113
根据式(9)计算：
Figure GDA00023400458300001114
联合上式可以建立以下成像模型：
Figure GDA00023400458300001115
Figure GDA00023400458300001116
和
Figure GDA00023400458300001117
分别为uj,k、vj,k关于参数u0,v0,f,θ1,θ2,θ3,
Figure GDA00023400458300001118
的函数；
4.4.2由步骤4.3给定星敏坐标系S参数矢量初值
Figure GDA00023400458300001119
畸变参数矢量初值
Figure GDA00023400458300001120
和星点在惯性坐标系i下的坐标
Figure GDA0002340045830000121
根据上式(10)给出的星点成像模型计算全部星点的图像坐标估计值，根据提取得到的星点图像坐标(uj,k,vj,k)计算星点图像坐标估计误差：
Figure GDA0002340045830000122
Figure GDA0002340045830000123
为利用式(10)计算得到的第j帧星图中第k个星点的图像坐标估计值，
Figure GDA0002340045830000124
为第j帧星图中第k个星点的图像坐标估计误差。
4.4.3估计星敏感器2的星敏坐标系S内外参数误差；
不考虑镜头畸变，星敏感器2的星敏坐标系S参数误差矢量
Figure GDA0002340045830000125
与星点图像坐标估计误差
Figure GDA0002340045830000126
满足如下关系：
Figure GDA0002340045830000127
其中，
Figure GDA0002340045830000128
为星敏坐标系S的参数矢量，
Figure GDA0002340045830000129
为X的估计误差矢量；
Figure GDA00023400458300001210
为
Figure GDA00023400458300001211
对参数矢量X的偏导数矢量；
Figure GDA00023400458300001212
为
Figure GDA00023400458300001213
对参数矢量X的偏导数矢量；
对步骤(4)中(4.2)提取的所有的星图的全部星点联立方程，完成步骤(4)中(4.3)；
Figure GDA00023400458300001214
其中：
Figure GDA00023400458300001215
N为全部星图的总帧数，M为第N帧星图的星点数；
将全部星点联立方程简写为Z＝H·ΔX，其中
Figure GDA0002340045830000131
采用最小二乘算法计算ΔX如下：
ΔX＝(HT·H)-1HT·Z (14)
4.4.4更新星敏感器2的星敏坐标系S内外参数，完成步骤(4)中(4.4)；
利用估计得到的星敏感器2内外参数误差更新星敏感器2的内外参数，更新方法如下：
Figure GDA0002340045830000132
然后将新的星敏感器2的星敏坐标系S参数矢量X赋值给内外参数初值
Figure GDA0002340045830000133
即
Figure GDA0002340045830000134
4.5星敏感器2的光学镜头畸变估计；
4.5.1给定星敏感器2的星敏坐标系S参数矢量
Figure GDA0002340045830000135
畸变参数矢量
Figure GDA0002340045830000136
和星点在惯性坐标系i下的坐标
Figure GDA0002340045830000137
根据式(10)给出的星点成像模型计算全部星点的图像坐标估计值
Figure GDA0002340045830000138
根据步骤(4.2)提取得到的星点图像坐标(uj,k,vj,k)计算星点图像坐标估计误差，完成步骤(4)中(4.5)：
Figure GDA0002340045830000139
4.5.2估计星敏坐标系S镜头畸变，优化内外参数；
星点图像坐标估计误差
Figure GDA00023400458300001310
满足如下关系：
Figure GDA00023400458300001311
其中：
Figure GDA00023400458300001312
为
Figure GDA00023400458300001313
对参数矢量X的偏导数矢量；
Figure GDA00023400458300001314
为
Figure GDA00023400458300001315
对参数矢量X的偏导数矢量，
Figure GDA00023400458300001316
ΔKd为畸变参数误差；
对于步骤4.2提取所有的星图提取的全部星点联立方程：
Figure GDA00023400458300001317
其中：
Figure GDA0002340045830000141
N为全部星图的总帧数，M为第N帧星图的星点数；
进一步将式(15)简写为Z＝H·ΔX+，其中
Figure GDA0002340045830000142
采用最小二乘算法计算ΔX如下：
ΔX＝(HT·H)-1HT·Z (19)
4.5.3更新星敏感器2星敏坐标系S内外参数和畸变参数；
利用估计得到的星敏感器2的星敏坐标系S内外参数误差更新星敏感器2的星敏坐标系S内外参数，更新方法如下：
Figure GDA0002340045830000143
然后将新的星敏感器2的星敏坐标系S参数矢量X赋值给内外参数初值
Figure GDA0002340045830000144
即
Figure GDA0002340045830000145
将新的星敏感器2的星敏坐标系S畸变参数矢量Kd赋值给畸变参数初值
Figure GDA0002340045830000146
即
Figure GDA0002340045830000147
4.6判断是否满足标定精度的要求；
利用标定得到的星敏感器2的星敏坐标系S参数矢量X、畸变参数矢量Kd和星点在惯性坐标系i下的坐标
Figure GDA0002340045830000148
根据式(10)给出的星点成像模型计算全部星点的图像坐标估计值
Figure GDA0002340045830000149
根据步骤4.2提取得到的星点图像坐标(uj,k,vj,k)计算星点图像坐标估计误差，完成步骤(4)中(4.6)：
Figure GDA00023400458300001410
如果
Figure GDA00023400458300001411
则跳转到步骤4.7，否则跳转到步骤5.4.2重新进行参数估计；其中max()为取最大值计算，||用于求矢量的模，ε为星点重构误差的阈值，ε的典型值为0.1像素；
4.7结束标定。

Claims (1)

1.一种多星图姿态关联的星敏感器内参数标定方法，采用以下步骤实现：
步骤(1)：标定开始；
(1.1)通过数据处理计算机控制陀螺组合体，使其开始进行姿态解算，并输出陀螺组合体相对于惯性坐标系的姿态，利用数据处理计算机采集陀螺组合体输出的姿态数据和每个数据所对应的UTC时间信息；
(1.2)利用数据处理计算机控制星敏感器使其开始拍摄恒星，并采集星敏感器拍摄得到的星图，记录每帧星图对应的UTC时间；
步骤(2)：标定数据采集；
在标定过程中通过手动或者二轴转台绕任意3个不同方向的转轴分别转动星惯组合系统，使星敏感器在不同的姿态下拍摄恒星；
步骤(3)：结束标定数据采集；控制陀螺组合体和星敏感器停止工作；
步骤(4)：标定数据处理；通过数据处理计算机运行算法进行标定数据处理；
所述步骤(4)中标定的数据处理具体的过程如下：
(4.1)标定数据的预处理；
(4.2)星点提取与全天时星图识别；
(4.3)给定星敏感器的初始参数矢量
Figure FDA0002364255860000012
和初始畸变参数
Figure FDA0002364255860000013
(4.4)估计星敏感器的主点、焦距以及外参数；
(4.5)星敏感器的光学镜头畸变估计；
(4.6)判断是否满足标定精度的要求；
(4.7)结束标定；
其特征在于：所述步骤(4)中(4.4)估计星敏感器的主点、焦距以及外参数和(4.5)星敏感器的光学镜头畸变估计算法如下：
(4.4.1)建立姿态帧关联星点成像模型；
建立星敏感器的星敏坐标系S星点成像模型如下：
Figure FDA0002364255860000011
其中，(uj,k,vj,k)为第j帧星图中第k个星点的图像坐标，(u0,v0)、f分别为待标定的星敏感器的主点、焦距，
Figure FDA0002364255860000021
为第j帧中的第k个星点在星敏坐标系S下的坐标，光学镜头畸变采用下式(6)给出的参数模型计算：
Figure FDA0002364255860000022
其中p1,p2,q1,q2,q3为待标定的星敏感器的镜头畸变参数，令畸变参数矢量为Kd＝[p1,p2,q1,q2,q3]T，rj,k为星点图像坐标(uj,k,vj,k)相对于主点(u0,v0)的距离；
对于第1帧星图，星点在星敏坐标系S与惯性坐标系i下的坐标存在如下映射关系：
Figure FDA0002364255860000023
其中
Figure FDA0002364255860000024
为第1帧星图中的第k个星点在惯性坐标系i下的坐标；
Figure FDA0002364255860000025
为第1帧星图时刻星敏坐标系S相对于惯性坐标系i的姿态矩阵，用欧拉角
Figure FDA0002364255860000026
表示如下：
Figure FDA0002364255860000027
第j帧的中的第k个星点在星敏坐标系S下的坐标通过第1帧星图的坐标递推得到：
Figure FDA0002364255860000028
其中，
Figure FDA0002364255860000029
为第j帧星图中的第k个星点在惯性坐标系下的坐标，
Figure FDA00023642558600000210
表示星敏坐标系S相对于陀螺组合体的安装矩阵，用欧拉角θ1,θ2,θ3表示如下：
Figure FDA00023642558600000211
Figure FDA00023642558600000212
为陀螺组合体坐标系b从t1时刻到tj时刻的姿态变换矩阵，由步骤(1.1)所述陀螺组合体解算得到t1和tj采样时刻陀螺组合体坐标系b相对于惯性坐标系i的姿态
Figure FDA00023642558600000213
和
Figure FDA00023642558600000214
根据式(9)计算：
Figure FDA00023642558600000215
联合上式建立以下成像模型：
Figure FDA0002364255860000031
Figure FDA0002364255860000032
和
Figure FDA0002364255860000033
分别为uj,k、vj,k关于参数
Figure FDA0002364255860000034
的函数；
(4.4.2)由给定星敏感器参数矢量初值
Figure FDA0002364255860000035
畸变参数矢量初值
Figure FDA0002364255860000036
和星点在惯性坐标系i下的坐标
Figure FDA0002364255860000037
根据上式(10)给出的星点成像模型计算全部星点的图像坐标估计值，根据提取得到的星点图像坐标(uj,k,vj,k)计算星点图像坐标估计误差：
Figure FDA0002364255860000038
Figure FDA0002364255860000039
为利用式(10)计算得到的第j帧星图中第k个星点的图像坐标估计值，
Figure FDA00023642558600000310
为第j帧星图中第k个星点的图像坐标估计误差；
(4.4.3)估计星敏感器内外参数误差：
不考虑镜头畸变，星敏感器的参数误差矢量
Figure FDA00023642558600000311
与星点图像坐标估计误差
Figure FDA00023642558600000312
满足如下关系：
Figure FDA00023642558600000313
其中，
Figure FDA00023642558600000314
为星敏感器的参数矢量，
Figure FDA00023642558600000315
为X的估计误差矢量；
Figure FDA00023642558600000316
为
Figure FDA00023642558600000317
对参数矢量X的偏导数矢量；
Figure FDA00023642558600000318
为
Figure FDA00023642558600000319
对参数矢量X的偏导数矢量；
对步骤(4)中(4.2)提取的所有的星图的全部星点联立方程；
Figure FDA0002364255860000041
其中：
Figure FDA0002364255860000042
为全部星图的总帧数，M为第N帧星图的星点数；
将全部星点联立方程简写为Z＝H·ΔX，其中
Figure FDA0002364255860000043
采用最小二乘算法计算ΔX如下：
ΔX＝(HT·H)-1HT·Z (14)
(4.4.4)更新星敏感器的内外参数；
利用估计得到的星敏感器内外参数误差更新星敏感器的内外参数，更新方法如下：
Figure FDA0002364255860000044
然后将新的星敏感器的参数矢量X赋值给内外参数初值
Figure FDA0002364255860000045
即
Figure FDA0002364255860000046
(4.5.1)给定星敏感器参数矢量
Figure FDA0002364255860000047
畸变参数矢量
Figure FDA0002364255860000048
和星点在惯性坐标系i下的坐标
Figure FDA0002364255860000049
根据式(10)给出的星点成像模型计算全部星点的图像坐标估计值
Figure FDA00023642558600000410
根据提取得到的星点图像坐标(uj,k,vj,k)计算星点图像坐标估计误差：
Figure FDA00023642558600000411
(4.5.2)估计星敏感器镜头畸变，优化内外参数；
星点图像坐标估计误差
Figure FDA00023642558600000412
满足如下关系：
Figure FDA0002364255860000051
其中，
Figure FDA0002364255860000052
为
Figure FDA0002364255860000053
对参数矢量X的偏导数矢量；
Figure FDA0002364255860000054
为
Figure FDA0002364255860000055
对参数矢量X的偏导数矢量，
Figure FDA0002364255860000056
ΔKd为畸变参数误差；
对于步骤(4.2)提取到的所有的星图的全部星点联立方程：
Figure FDA0002364255860000057
其中：
Figure FDA0002364255860000058
N为全部星图的总帧数，M为第N帧星图的星点数；
将上式全部星点联立方程简写为Z＝H·ΔX+，其中
Figure FDA0002364255860000059
采用最小二乘算法计算ΔX如下：
ΔX＝(HT·H)-1HT·Z (19)
(4.5.3)更新星敏感器内外参数和畸变参数；
利用估计得到的星敏感器内外参数误差更新星敏感器的内外参数，更新方法如下：
Figure FDA0002364255860000061
然后将新的星敏感器的参数矢量X赋值给内外参数初值
Figure FDA0002364255860000062
即
Figure FDA0002364255860000063
将新的星敏感器的畸变参数矢量Kd赋值给畸变参数初值
Figure FDA0002364255860000064
所述步骤(4.6)中判断是否满足标定精度的要求：利用标定得到的星敏感器参数矢量X、畸变参数矢量Kd和星点在惯性坐标系i下的坐标
Figure FDA0002364255860000065
根据式(10)给出的星点成像模型计算全部星点的图像坐标估计值
Figure FDA0002364255860000066
根据提取得到的星点图像坐标(uj,k,vj,k)计算星点图像坐标估计误差：
Figure FDA0002364255860000067
如果
Figure FDA0002364255860000068
则跳转到步骤(4.7)，否则跳转到步骤(4.4.2)重新进行参数估计，其中max()为取最大值函数，| |用于求矢量的模，ε为星点重构误差的阈值。
CN201710581215.0A 2017-07-17 2017-07-17 一种多星图姿态关联的星敏感器内参数标定方法 Active CN107449444B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710581215.0A CN107449444B (zh) 2017-07-17 2017-07-17 一种多星图姿态关联的星敏感器内参数标定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710581215.0A CN107449444B (zh) 2017-07-17 2017-07-17 一种多星图姿态关联的星敏感器内参数标定方法

Publications (2)

Publication Number Publication Date
CN107449444A CN107449444A (zh) 2017-12-08
CN107449444B true CN107449444B (zh) 2020-04-10

Family

ID=60487725

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710581215.0A Active CN107449444B (zh) 2017-07-17 2017-07-17 一种多星图姿态关联的星敏感器内参数标定方法

Country Status (1)

Country Link
CN (1) CN107449444B (zh)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107588785B (zh) * 2017-09-12 2019-11-05 中国人民解放军国防科技大学 一种考虑像点误差的星敏感器内外参数简化标定方法
CN108592945B (zh) * 2018-03-27 2020-08-21 中国人民解放军国防科技大学 一种惯性/天文组合系统误差的在线标定方法
CN108645401B (zh) * 2018-04-03 2020-05-22 中国人民解放军国防科技大学 基于姿态关联图像叠加的全天时星敏感器星点提取方法
CN109506656B (zh) * 2018-11-28 2020-11-10 上海航天控制技术研究所 一种高精度在轨姿态信息下传还原方法
CN109459065B (zh) * 2018-12-26 2020-06-19 长光卫星技术有限公司 一种基于卫星惯性空间旋转姿态的陀螺安装矩阵标定方法
CN109798921B (zh) * 2019-02-22 2022-07-19 中国科学院光电技术研究所 一种星敏感器内方位元素室内校准方法
CN111402300B (zh) * 2020-04-21 2022-09-20 中国科学院光电技术研究所 一种基于双谱域主成分分析的高动态星敏感器运动参数估计方法
CN112254743B (zh) * 2020-10-15 2024-05-31 长春工业大学 一种基于星角距相减的星敏感器在轨标定方法
CN112882483B (zh) * 2021-01-12 2022-03-04 北京控制工程研究所 星敏感器在轨标定方法、装置和存储介质
CN112989637B (zh) * 2021-05-06 2021-07-23 中国人民解放军国防科技大学 基于分布折算、近似和综合的星敏系统可靠性评估方法
CN113514055B (zh) * 2021-07-09 2022-10-28 北京航空航天大学 一种地基星敏感器大气折射与对地姿态的联合估计方法
CN114663492B (zh) * 2022-02-16 2024-02-09 西北工业大学 一种利用事件相机作为星敏感器的飞行器姿态确定方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101013033A (zh) * 2006-03-21 2007-08-08 北京航空航天大学 一种基于无偏差带的星敏感器地面校准方法
US8019544B2 (en) * 2005-01-03 2011-09-13 The Boeing Company Real-time refinement method of spacecraft star tracker alignment estimates
CN105659815B (zh) * 2007-08-08 2012-01-04 北京航空航天大学 一种星敏感器动态校准装置和校准方法
CN103674023A (zh) * 2013-12-26 2014-03-26 中国人民解放军国防科学技术大学 一种基于陀螺精确角度关联的星敏感器动态测姿方法
CN105737858A (zh) * 2016-05-04 2016-07-06 北京航空航天大学 一种机载惯导系统姿态参数校准方法与装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8019544B2 (en) * 2005-01-03 2011-09-13 The Boeing Company Real-time refinement method of spacecraft star tracker alignment estimates
CN101013033A (zh) * 2006-03-21 2007-08-08 北京航空航天大学 一种基于无偏差带的星敏感器地面校准方法
CN105659815B (zh) * 2007-08-08 2012-01-04 北京航空航天大学 一种星敏感器动态校准装置和校准方法
CN103674023A (zh) * 2013-12-26 2014-03-26 中国人民解放军国防科学技术大学 一种基于陀螺精确角度关联的星敏感器动态测姿方法
CN105737858A (zh) * 2016-05-04 2016-07-06 北京航空航天大学 一种机载惯导系统姿态参数校准方法与装置

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
Star sensor calibration based on integrated modelling with intrinsic and extrinsic parameters;Xinguo Wei 等;《Measurement》;20140930;第55卷;第117-125页 *
星敏感器外场观星标定及检验方法研究;姜文英 等;《计量学报》;20160531;第37卷(第3期);第251-254页 *
星敏感器模型参数分析与校准方法研究;郝雪涛 等;《光电工程》;20050331;第32卷(第3期);第5-8页 *
高精度全天时星敏感器关键技术研究;何家维;《中国博士学位论文全文数据库 工程科技Ⅱ辑》;20140215(第02期);第C031-6页 *
高精度星敏感器的标定;乔培玉 等;《红外与激光工程》;20121031;第41卷(第10期);第2779-2784页 *

Also Published As

Publication number Publication date
CN107449444A (zh) 2017-12-08

Similar Documents

Publication Publication Date Title
CN107449444B (zh) 一种多星图姿态关联的星敏感器内参数标定方法
CN111947652B (zh) 一种适用于月球着陆器的惯性/视觉/天文/激光测距组合导航方法
CN110009681B (zh) 一种基于imu辅助的单目视觉里程计位姿处理方法
CN106780699B (zh) 一种基于sins/gps和里程计辅助的视觉slam方法
CN109029433B (zh) 一种移动平台上基于视觉和惯导融合slam的标定外参和时序的方法
CN108592945B (zh) 一种惯性/天文组合系统误差的在线标定方法
Veth Fusion of imaging and inertial sensors for navigation
CN101598556B (zh) 一种未知环境下无人机视觉/惯性组合导航方法
CN102741706B (zh) 地理参照图像区域的方法
CN107607127B (zh) 一种基于外场的星敏感器内部参数标定及精度快速验证系统
CN107796391A (zh) 一种捷联惯性导航系统/视觉里程计组合导航方法
CN109459059B (zh) 一种星敏感器外场转换基准测定系统及方法
CN109520476B (zh) 基于惯性测量单元的后方交会动态位姿测量系统及方法
CN113551668B (zh) 一种航天器惯性/恒星星光矢量/星光折射组合导航方法
CN111665512A (zh) 基于3d激光雷达和惯性测量单元的融合的测距和绘图
CN111307146A (zh) 一种基于双目相机和imu的虚拟现实头戴显示设备定位系统
CN112797985A (zh) 基于加权扩展卡尔曼滤波的室内定位方法及室内定位系统
CN108444468B (zh) 一种融合下视视觉与惯导信息的定向罗盘
Lo et al. The direct georeferencing application and performance analysis of uav helicopter in gcp-free area
CN114964276A (zh) 一种融合惯导的动态视觉slam方法
CN113218577A (zh) 一种星敏感器星点质心位置精度的外场测量方法
CN117053797A (zh) 基于多目视觉下的大气偏振导航方法
CN114001756B (zh) 一种小视场星敏感器外场地面寻星方法
CN109064510B (zh) 一种全站仪及其恒星图像的星点质心提取方法
CN112197765B (zh) 一种实现水下机器人精细导航的方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant