CN110006462B - 基于奇异值分解的星敏感器在轨标定方法 - Google Patents

基于奇异值分解的星敏感器在轨标定方法 Download PDF

Info

Publication number
CN110006462B
CN110006462B CN201910432037.4A CN201910432037A CN110006462B CN 110006462 B CN110006462 B CN 110006462B CN 201910432037 A CN201910432037 A CN 201910432037A CN 110006462 B CN110006462 B CN 110006462B
Authority
CN
China
Prior art keywords
star
model
vector
matrix
singular
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
CN201910432037.4A
Other languages
English (en)
Other versions
CN110006462A (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.)
Changchun University of Technology
Original Assignee
Changchun University of 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 Changchun University of Technology filed Critical Changchun University of Technology
Priority to CN201910432037.4A priority Critical patent/CN110006462B/zh
Publication of CN110006462A publication Critical patent/CN110006462A/zh
Application granted granted Critical
Publication of CN110006462B publication Critical patent/CN110006462B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C25/00Manufacturing, 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

本发明涉及一种基于奇异值分解的星敏感器在轨标定方法,属于天文导航领域。包括如下步骤:利用坐标变换过程中奇异值不变的特性建立校准模型;利用可观测性对模型进行优化;最优星组合模型的选择;奇异值的选择;基于扩展卡尔曼滤波器EKF估计相机参数。与传统的角距离(AD)方法相比,本发明方法的精度和收敛速度好,SVD方法的每帧平均标定时间明显短于AD方法。这意味着SVD方法的计算成本比传统方法明显缩小,这对于资源有限的星敏感器具有非凡的意义。

Description

基于奇异值分解的星敏感器在轨标定方法
技术领域
本发明涉及天文导航领域,特别涉及一种标定方法,尤指一种基于奇异值分解的星敏感器在轨标定方法。
背景技术
星敏感器是一种导航系统,通过对恒星的观测,获取载体的姿态信息。它是目前最精确的光学姿态传感器。由于其导航精度高、自主性强、无累积误差,受到航空航天工业的青睐。作为航天器的“眼睛”,星敏感器的精度直接决定了航天器的性能。然而,星敏感器是一种光学器件,其精度取决于成像质量和光学参数(包括焦距、主点和畸变)的精度。因此,标定是星敏感器的关键技术之一。
星敏感器的标定方法可分为地面标定和在轨标定两类。一般地,地面标定方法需要积累大量的测量数据,因此该方法往往依赖于一个固定平台和复杂的实验机制,标定成本高。此外,星敏感器的实际工作环境与校准环境有着明显的不同,空间参数可能发生变化。在大多数情况下,基于运行中的观测数据进行在轨标定是必要的,以便实时更新参数,保持星敏感器的精度。目前星敏感器在轨标定方法主要依赖于星与星之间的角距离,很少有研究致力于寻找新的标定参考,这种方法数据存储量较大、计算效率较低。
发明内容
本发明的目的在于提供一种基于奇异值分解的星敏感器在轨标定方法,解决了现有技术存在的上述问题。针对传统方法的不足,本发明提供一种基于奇异值分解的星敏感器在轨标定方法,在星敏感器标定中具有良好的应用效果,尤其适用于计算资源有限的星敏感器。本发明采用不变奇异值做为定标参考,基于星矢量在坐标变换下奇异值不变的在轨标定方法,与传统的基于角距离的方法相比,本发明的计算成本明显降低。
本发明的上述目的通过以下技术方案实现:
基于奇异值分解的星敏感器在轨标定方法,包括如下步骤:
步骤1)、利用坐标变换过程中奇异值不变的特性建立校准模型;
步骤2)、利用可观测性对模型进行优化;
步骤2.1)、最优星组合模型的选择;
步骤2.2)、奇异值的选择;
步骤3)、基于扩展卡尔曼滤波器EKF估计相机参数。
步骤1)所述的利用坐标变换过程中奇异值不变的特性建立校准模型是:
通过相机模型,在图像坐标系中给定观测到的恒星质心的位置pd=[ud,vd]T,求出相机坐标系w=[x,y,z]T中对应的恒星矢量,并将其关系表示为:
w=F(s,f,u0,v0,k1,k2,pd) (1)
其中F(·)是带有失真的反投影函数,s是纵横比,f是焦距,[u0,v0]T为图像坐标系中主点的坐标,k1、k2为畸变系数;
定义wi是惯性坐标系中的一个导航星单位矢量,vi是相机坐标系中的一个观测星矢量;这两个坐标系之间的变换是:
W=CV (2)
其中W和V是列向量矩阵
W=[w1 w2……wN]3×N (3)
V=[v1 v2……vN]3×N (4)
C为姿态矩阵,表示从惯性坐标系到相机坐标系的变换,因此C为正交矩阵;
利用奇异值分解法,将矩阵W和V分解为
Figure SMS_1
Figure SMS_2
Pw和Pv是左奇异向量pwi和pvi的3×3正交矩阵(i=1,2,3),Qv和Qw是右奇异向量qvi和qwi的N×N正交矩阵(i=1,2,3),∑v和∑w是3×N对角矩阵,对角元素是V和W的奇异值σvi和σwi(i=1,2,3);对于视场中三个以上的不同恒星,有3个非零奇异值,SVD是唯一的;
将公式(2)乘以WT,可以得到
WWT=CVVTCT (7)
将式(5)和(6)带入式(7)得到
Figure SMS_3
以及
Figure SMS_4
其中Sw和Sv是具有特征值
Figure SMS_5
Figure SMS_6
的WWT和VVT的对角矩阵;由于WWT和VVT是正定对称矩阵,C是正交矩阵,式(7)是相似变换;因此,WWT和VVT的特征值相等,即
Figure SMS_7
因此,W和V的正奇异值相等:
σwi=σvi,i=1,2,3 (11)
假设SV(·)是奇异值求解算子,公式(11)可以表示为
σvi=SV(V)=σwi=SV(W),i=1,2,3 (12)
将式(1)带入(12)得到
σvi=σwi=SV(V)=SV(F(s,f,u0,v0,k1,k2,Pd)),i=1,2,3 (13)
其中Pd=[pd1 pd2……pdN]3×N是图像坐标系中观察到的星坐标的集合;
恒星识别后,观测到的恒星坐标Pd与星表中相应的恒星矢量V相互匹配;因此,根据公式(13),σvi可以由星向量V得到,也可以由摄像机参数和观测到的星坐标Pd计算;星表的精度很高,用V求出的观测量σvi具有很好的精度。
步骤2)所述的利用可观测性对模型进行优化是:
可观测性是评价系统可行性的指标,即在不同的模型下,相同的输入推导可能导致不同的输出推导;如果输出推导的幅度较大,则可观测性更好,系统更可行,反之亦然;根据可观测性的定义,可以得到
δzk=HkWδx (13)
其中δx是所有参数具有相同推导的输入推导向量,δzk是输出推导,Hk是雅可比矩阵;由于实际中不同参数的精度差异很大,这意味着δx不能代表不同参数的实际推导;δx应根据不同精度的大小进行加权,W是对角加权矩阵,其元素是根据现场实验中的星敏感器得出的典型参数精度;
所以可观测性矩阵是
H'k=HkW (14)
再次使用奇异值分解,利用可观测矩阵的奇异值分解,公式(13)可表示为
δzk=PkkQkδx (15)
其中Pk和Qk分别是左奇异向量和右奇异向量的正交矩阵;为了确保H’k是可观测的,将∑k定义为6×N对角矩阵,对角元素是非零奇异值σi(i=1~6);
由于Pk和Qk是正交矩阵,可计算
Figure SMS_8
其中||δxk||2和||δzk||2分别为δx和δzk的2范数;
输出推导的下确界为
Figure SMS_9
其中σmin是可观测矩阵的最小奇异值(MSV);很明显,σmin越大,输出推导的最小值越大,可观测性越好;因此,采用σmin作为评价系统性能的指标,寻找合适的标定模型。
步骤2.1)所述的最优星组合模型的选择是:
对于超过三颗星,有三个非零奇异值,因此组合可以由三颗星、四颗星等等组成;如果所有的组合都被用来构成观测模型,计算量是巨大的;此外,组合不独立,校准信息不随组合数量的增加而增加;
设视场中有N颗星星,考虑如下几种模型:
模型1:三颗星组成一个组合,如1-2-3、2-3-4、3-4-5;
模型2:
Figure SMS_10
个星,
Figure SMS_11
为向下取整,构成一个组合;
模型3:N-1个星构成一个组合;
模型4:与其他模型不同,该模型结合了不同数量的星形组合,如1-2-3、1-2-3-4、1-2-3-4-5;
根据可观测性分析,模型4的最小奇异值总是最大的,这意味着模型4的可观测性比其他模型好,所以将模型4作为最佳组合模型。
步骤2.2)所述的奇异值的选择是:
三个奇异值对状态变量变化的敏感度是完全不同的;σ12和σ3是三个具有降序的奇异值,σ1的可观测性比σ23,差,因此我们应该考虑这三个奇异值是否都适用于摄像机参数的估计。
矩阵W的奇异值分解的物理解释如图1所示。黑色矢量表示恒星矢量wi;根据奇异值分解的定义,左奇异向量pwi(i=1,2,3)是单位向量,因此W的最大奇异值的平方为
Figure SMS_12
其中pmax表示与最大奇异值σ1相关的左奇异向量;
公式(18)表明,pmax是使每个恒星矢量wi的投影之和达到最大值的矢量;同样,与最小奇异值σ3相关联的pmin是使每个星向量wi投影到pmin上之和最小的向量;然后,与中间奇异值相关联的pmiddle垂直于pmax和pmin生成的平面;
由于pmax和wi是单位向量,因此公式(18)可以重写为
Figure SMS_13
其中αi是wi与pmax的夹角,它也给出了星矢量wi对pmax的投影;因此,σ1的推导在sinαi的大小上,由于pmiddle和pmin生成的平面垂直于pmax,因此σ2和σ3与sinαi成比例,σ2和σ3的推导在cosαi的大小上,对于较小的视场,几乎所有sinαi都小于cosαi,因此σ1的推导小于σ2和σ3,这与可观测性分析是一致的;因此,为了减少计算量,采用σ2和σ3作为观测量。
步骤3)所述的基于扩展卡尔曼滤波器EKF估计相机参数,具体是:
状态转换和测量模型是
xk=I6×6·xk-1 (20)
zk=h(xk)+nc,i=1,2,3 (21)
其中xk=[s,f,u0,v0,k1,k2]T,xk和xk-1分别是标记为k和k-1的恒星图像的状态,相机参数为常数,I6×6是一个单位矩阵,zk=[σ23]T可以用导航星矢量V计算,h(xk)是 SV(F(s,f,u0,v0,k1,k2,Pd))的简化表示,nc是由噪声引起的测量误差;
从地面定标得到的恒星协方差P0和参数x0的初始估计开始,逐帧处理校准;对于第k 个星图像,EKF预测方程为
Figure SMS_14
Figure SMS_15
其中Q是先验估计误差的协方差矩阵;
EKF更新方程为
Figure SMS_16
Figure SMS_17
Figure SMS_18
其中R是观测噪声的协方差矩阵;Hk是观测的雅可比矩阵
Figure SMS_19
由于h模型是复杂的,所以采用数值微分法来计算雅可比矩阵。
本发明的有益效果在于:与传统的角距离(AD)方法相比,奇异值分解SVD方法的精度和收敛速度好,SVD方法的每帧平均标定时间明显短于AD方法。这意味着SVD方法的计算成本比传统方法明显缩小,这对于资源有限的星敏感器具有非凡的意义。
附图说明
此处所说明的附图用来提供对本发明的进一步理解,构成本申请的一部分,本发明的示意性实例及其说明用于解释本发明,并不构成对本发明的不当限定。
图1为矩阵w的奇异值分解物理解释;
图2为本发明的星组合模型4的选取示意图。
具体实施方式
下面结合附图进一步说明本发明的详细内容及其具体实施方式。
参见图1及图2所示,本发明的基于奇异值分解的星敏感器在轨标定方法,具体步骤如下:
1.建立基于奇异值的校准模型:
利用相机模型,在图像坐标系中给定观测到的恒星质心的位置pd=[ud,vd]T,可以求出相机坐标系w=[x,y,z]T中对应的恒星矢量,并将其关系表示为
w=F(s,f,u0,v0,k1,k2,pd) (1)
其中F(·)是带有失真的反投影函数,s是纵横比,f是焦距,[u0,v0]T为图像坐标系中主点的坐标,k1、k2为畸变系数。
定义wi是惯性坐标系中的一个导航星单位矢量,vi是相机坐标系中的一个观测星矢量。这两个坐标系之间的变换是
W=CV (2)
其中W和V是列向量矩阵,用模型4(图2)选择星组合模型。
W=[w1 w2 wN]3×N (3)
V=[v1 v2 vN]3×N (4)
C为姿态矩阵,表示从惯性坐标系到相机坐标系的变换,因此C为正交矩阵。
利用奇异值分解法,将矩阵W和V分解为
Figure SMS_20
Figure SMS_21
Pw和Pv是左奇异向量pwi和pvi的3×3正交矩阵(i=1,2,3),Qv和Qw是右奇异向量qvi和qwi的N×N正交矩阵(i=1,2,3),∑v和∑w是3×N对角矩阵,对角元素是V和W的奇异值σvi和σwi(i=1,2,3)。对于视场中三个以上的不同恒星,有3个非零奇异值,SVD是唯一的。奇异值具有以下特性。
将公式(1)乘以WT,可以得到
WWT=CVVTCT (7)
将式(5)和(6)带入式(7)得到
Figure SMS_22
以及
Figure SMS_23
其中Sw和Sv是具有特征值
Figure SMS_24
Figure SMS_25
的WWT和VVT的对角矩阵。由于WWT和 VVT是正定对称矩阵,C是正交矩阵,式(6)是相似变换。因此,WWT和VVT的特征值相等,即
Figure SMS_26
因此,W和V的正奇异值相等:
σwi=σvi,i=1,2,3 (11)
假设SV(·)是奇异值求解算子,公式(11)可以表示为
σvi=SV(V)=σwi=SV(W),i=1,2,3 (12)
将式(1)带入(12)得到
σvi=σwi=SV(V)=SV(F(s,f,u0,v0,k1,k2,Pd)),i=1,2,3 (13)
其中Pd=[pd1 pd2……pdN]3×N是图像坐标系中观察到的星坐标的集合。
恒星识别后,观测到的恒星坐标Pd与星表中相应的恒星矢量V相互匹配。因此,根据公式(13),σvi可以由星向量V得到,也可以由摄像机参数和观测到的星坐标Pd计算。星表的精度很高,用V求出的观测量σvi具有很好的精度。
2.利用上述模型进行在轨标定
在某一时刻下,由光学视场捕获恒星目标,并成像在探测器上,通过对探测器上的图像进行星图预处理及质心计算获取恒星在图像中的质心位置,即标定过程中的恒星图像坐标位置。利用图像中恒星的位置,结合地面初始标定获取的光学系统参数信息,可以得到相机坐标系中星点坐标的粗略位置,通过星图识别与星表中的恒星矢量位置比对,即可获取与图像中恒星位置相对应的空间坐标位置。此时已知恒星的图像坐标和空间坐标,即可进行在轨标定工作。
基于扩展卡尔曼滤波器(EKF)估计相机参数。状态转换和测量模型是
xk=I6×6·xk-1 (14)
zk=h(xk)+nc,i=1,2,3 (15)
其中xk=[s,f,u0,v0,k1,k2]T,xk和xk-1分别是标记为k和k-1的恒星图像的状态,假设相机参数为常数,I6×6是一个单位矩阵,zk=[σ23]T可以用导航星矢量V计算,h(xk)是SV(F(s,f,u0,v0,k1,k2,Pd))的简化表示,nc是由噪声引起的测量误差。
从地面标定得到的噪声协方差P0和参数x0的初始估计开始,我们逐帧处理校准。对于第k个星图像,EKF预测方程为
Figure SMS_27
Figure SMS_28
其中Q是先验估计误差的协方差矩阵。
EKF更新方程为
Figure SMS_29
Figure SMS_30
Figure SMS_31
其中R是观测噪声的协方差矩阵。Hk是观测的雅可比矩阵
Figure SMS_32
由于h模型是复杂的,所以我们采用数值微分法来计算雅可比矩阵。
3.仿真实验与分析
采用19.14°×11.18°视场、1920×1080像素阵列星敏感器,以2Hz更新率进行模拟。其他模拟参数如上表典型参数精度所示。模拟数据由三组数据组成:惯性系S3D中的3D星矢量、图像系S2D中相应的2D星坐标和具有正态分布噪声的2D星坐标集(标准差为0.5像素),即S’2D,对2500幅图像的S3D、S2D和S’2D进行了仿真,这些图像分为两组,第一组 2400幅图像用于标定(根据实验,我们发现2400幅图像可以保证标定的收敛性),最后100 幅图像用于标定性能的评估。
为了充分评价标定方法的性能,提出了以下评价标准。由于奇异值的残差不能直接反映标定结果对星敏感器性能的影响,本发明采用了星间角距离残差来评价在轨标定的性能。然后,采用两个准则来评价校准的性能。
准则A:该准则基于无误差的仿真数据。通过对摄像机参数的估计,可以将S2D反投影到摄像机坐标系上,并计算出它们之间的角度距离。利用这些角距离和S3D计算出的相应角距离得到了残差,然后计算出了每幅恒星图像中角距离的均方根误差(RMSE)(A|δRMSE)。当估计收敛时,最终的评价指标是最近100幅图像中A|δRMSE的均值
Figure SMS_33
和标准差
Figure SMS_34
准则B:该准则与准则A的区别在于,该模型使用有噪声的2D模拟数据S’2D来计算每个恒星图像中角距离(B|δRMSE)的RMSE。最后100幅图像以B|δRMSE的平均值
Figure SMS_35
和标准差
Figure SMS_36
表示标准B的评价指标。
首先对星组合模型的性能进行测试,恒星星表是通过从Tycho-2目录中选取5.5级的恒星形成的。初始值设置为:s=1,f=15.5mm,u0=960,v0=540,k1=0,k2=0。校准结果如下表。
Figure SMS_37
由于参数的尺度不同,上表中所示参数的推导表示为分数变化。结果表明,模型2和模型3得到的参数s,u0,v0,k2的精度优于模型1;模型1和模型2得到的f的精度优于模型3;模型1得到的k1的精度优于模型2和模型3。基于这些结果,我们得出结论:对于不同的参数,模型1、模型2和模型3各有优缺点。然而,模型4标定的参数精度始终处于中间水平,这意味着模型4结合了其他组合模型的特性。因此我们选择模型4进行其他实验。
接下来测试不同奇异值组合的性能,利用奇异值的所有组合构造雅可比矩阵,校正结果如下表。
Figure SMS_38
根据上表的结果,三种最佳组合是{σ123},{σ23}和{σ2}。即σ2和σ3的可观测性优于σ1。特别是,与σ2相关的雅可比矩阵的可观测性优于σ3,因此仅用σ2获得的标定结果是足够好的。但通过多次实验,我们发现{σ2}的性能不稳定,在本实验中,相对较高的标准差
Figure SMS_39
给出了稳定性问题的一些指示。这可能是由于σ1和σ3中的一些信息有助于校准,因此建议使用{σ123},{σ23}。{σ123}和{σ23}的校准精度相似,但{σ23}的经过时间短于{σ1, σ23},这意味着{σ23}的计算小于{σ123},因此在以下试验中采用{σ23}的组合。
建立三个实验组,实验组之间的差异是恒星传感器的极限视觉星等,从而导致视场中的恒星数目不同。各组中的平均恒星数分别为31.4、19.1和7.7,极限视星等分别为6、5.5和 4.6。结果见下表。
Figure SMS_40
在标准A下,SVD方法的性能接近甚至优于AD方法,而且SVD方法更稳定。然而,在标准B下,SVD方法的性能总是比AD方法差。这是因为AD方法以角距离为优化估计对象,因此综合考虑摄像机参数和中心噪声的标定效果较好。必须指出的是,星敏感器的姿态估计是基于星向量而不是角距离,因此摄像机参数的精度更为重要。因此,SVD方法的精度与AD方法相当。两种方法的收敛速度也很相似,这是针对第2组的。然而,SVD方法的运行时间比AD方法短,特别是对于视场中恒星较多的情况。对于第一组,SVD方法的计算速度比AD方法提高了94.48%。
以上所述仅为本发明的优选实例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡对本发明所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于奇异值分解的星敏感器在轨标定方法,其特征在于:包括如下步骤:
步骤1)、利用坐标变换过程中奇异值不变的特性建立校准模型;
步骤2)、利用可观测性对模型进行优化;
步骤2.1)、最优星组合模型的选择;
步骤2.2)、奇异值的选择;
步骤3)、基于扩展卡尔曼滤波器EKF估计相机参数;
其中,步骤1)所述的利用坐标变换过程中奇异值不变的特性建立校准模型是:
1.1)、建立基于奇异值的校准模型:
通过相机模型,在图像坐标系中给定观测到的恒星质心的位置pd=[ud,vd]T,求出相机坐标系w=[x,y,z]T中对应的恒星矢量,并将其关系表示为:
w=F(s,f,u0,v0,k1,k2,pd) (1)
其中F(·)是带有失真的反投影函数,s是纵横比,f是焦距,[u0,v0]T为图像坐标系中主点的坐标,k1、k2为畸变系数;
定义wi是惯性坐标系中的一个导航星单位矢量,vi是相机坐标系中的一个观测星矢量;这两个坐标系之间的变换是:
W=CV (2)
其中W和V是列向量矩阵
W=[w1 w2 … wN]3×N (3)
V=[v1 v2 … vN]3×N (4)
C为姿态矩阵,表示从惯性坐标系到相机坐标系的变换,因此C为正交矩阵;
利用奇异值分解法,将矩阵W和V分解为
Figure FDA0004053802380000011
Figure FDA0004053802380000021
Pw和Pv是左奇异向量pvi和pwi的3×3正交矩阵(i=1,2,3),Qv和Qw是右奇异向量qvi和qwi的N×N正交矩阵(i=1,2,3),∑v和∑w是3×N对角矩阵,对角元素是V和W的奇异值σvi和σwi(i=1,2,3);对于视场中三个以上的不同恒星,有3个非零奇异值,SVD是唯一的;
将公式(2)乘以WT,可以得到
WWT=CVVTCT (7)
将式(5)和(6)带入式(7)得到
Figure FDA0004053802380000022
以及
Figure FDA0004053802380000023
其中
Figure FDA0004053802380000024
Figure FDA0004053802380000025
是具有特征值
Figure FDA0004053802380000026
Figure FDA0004053802380000027
的VVT和WWT的对角矩阵;由于VVT和WWT是正定对称矩阵,C是正交矩阵,式(7)是相似变换;因此,VVT和WWT的特征值相等,即
Figure FDA0004053802380000028
因此,W和V的正奇异值相等:
σwi=σvi,i=1,2,3 (11)
假设SV(·)是奇异值求解算子,公式(11)可以表示为
σvi=SV(V)=σwi=SV(W),i=1,2,3 (12)
将式(1)带入(12)得到
σvi=σwi=SV(V)=SV(F(s,f,u0,v0,k1,k2,Pd)),i=1,2,3 (13)
其中Pd=[pd1 pd2 …… pdN]3×N是图像坐标系中观察到的星坐标的集合;
恒星识别后,观测到的恒星坐标Pd与星表中相应的恒星矢量V相互匹配;因此,根据公式(13),σvi可以由星向量V得到,也可以由摄像机参数和观测到的星坐标Pd计算;星表的精度很高,用V求出的观测量σvi具有很好的精度;
1.2)、利用上述模型进行在轨标定
在某一时刻下,由光学视场捕获恒星目标,并成像在探测器上,通过对探测器上的图像进行星图预处理及质心计算获取恒星在图像中的质心位置,即标定过程中的恒星图像坐标位置;利用图像中恒星的位置,结合地面初始标定获取的光学系统参数信息,可以得到相机坐标系中星点坐标的粗略位置,通过星图识别与星表中的恒星矢量位置比对,即可获取与图像中恒星位置相对应的空间坐标位置;此时已知恒星的图像坐标和空间坐标,即可进行在轨标定工作;
基于扩展卡尔曼滤波器(EKF)估计相机参数;状态转换和测量模型是
xk=I6×6·xk-1 (14)
zk=h(xk)+nc,i=1,2,3 (15)
其中xk=[s,f,u0,v0,k1,k2]T,xk和xk-1分别是标记为k和k-1的恒星图像的状态,假设相机参数为常数,I6×6是一个单位矩阵,zk=[σ23]T可以用导航星矢量V计算,h(xk)是SV(F(s,f,u0,v0,k1,k2,Pd))的简化表示,nc是由噪声引起的测量误差;从地面标定得到的噪声协方差P0和参数x0的初始估计开始,逐帧处理校准;对于第k个星图像,EKF预测方程为
Figure FDA0004053802380000031
Figure FDA0004053802380000032
其中Q是先验估计误差的协方差矩阵;
EKF更新方程为
Figure FDA0004053802380000033
Figure FDA0004053802380000034
Figure FDA0004053802380000035
其中R是观测噪声的协方差矩阵;Hk是观测的雅可比矩阵
Figure FDA0004053802380000041
由于h模型是复杂的,所以我们采用数值微分法来计算雅可比矩阵;
步骤2)所述的利用可观测性对模型进行优化是:
可观测性是评价系统可行性的指标,即在不同的模型下,相同的输入推导可能导致不同的输出推导;如果输出推导的幅度较大,则可观测性更好,系统更可行,反之亦然;根据可观测性的定义,可以得到
δzk=HkWδx (22)
其中δx是所有参数具有相同推导的输入推导向量,δzk是输出推导,Hk是雅可比矩阵;由于实际中不同参数的精度差异很大,这意味着δx不能代表不同参数的实际推导;δx应根据不同精度的大小进行加权,W是对角加权矩阵,其元素是根据现场实验中的星敏感器得出的典型参数精度;
所以可观测性矩阵是
H'k=HkW (23)
再次使用奇异值分解,利用可观测矩阵的奇异值分解,公式(22)可表示为
δzk=PkkQkδx (24)
其中Pk和Qk分别是左奇异向量和右奇异向量的正交矩阵;为了确保H’k是可观测的,将∑k定义为6×N对角矩阵,对角元素是非零奇异值σi(i=1~6);
由于Pk和Qk是正交矩阵,可计算
Figure FDA0004053802380000042
其中||δxk||2和||δzk||2分别为δx和δzk的2范数;
输出推导的下确界为
Figure FDA0004053802380000043
其中σmin是可观测矩阵的最小奇异值(MSV);很明显,σmin越大,输出推导的最小值越大,可观测性越好;因此,采用σmin作为评价系统性能的指标,寻找合适的标定模型;
步骤2.2)所述的奇异值的选择是:
三个奇异值对状态变量变化的敏感度是完全不同的;σ12和σ3是三个具有降序的奇异值,σ1的可观测性比σ23差;
矩阵W的奇异值分解的物理解释,黑色矢量表示恒星矢量wi;根据奇异值分解的定义,左奇异向量pwi(i=1,2,3)是单位向量,因此W的最大奇异值的平方为
Figure FDA0004053802380000051
其中pmax表示与最大奇异值σ1相关的左奇异向量;
公式(27)表明,pmax是使每个恒星矢量wi的投影之和达到最大值的矢量;同样,与最小奇异值σ3相关联的pmin是使每个星向量wi投影到pmin上之和最小的向量;然后,与中间奇异值相关联的pmiddle垂直于pmax和pmin生成的平面;
由于pmax和wi是单位向量,因此公式(27)可以重写为
Figure FDA0004053802380000052
其中αi是wi与pmax的夹角,它也给出了星矢量wi对pmax的投影;因此,σ1的推导在sinαi的大小上,由于pmiddle和pmin生成的平面垂直于pmax,因此σ2和σ3与sinαi成比例,σ2和σ3的推导在cosαi的大小上,对于较小的视场,所有sinαi都小于cosαi,因此σ1的推导小于σ2和σ3,这与可观测性分析是一致的;因此,为了减少计算量,采用σ2和σ3作为观测量。
2.根据权利要求1所述的基于奇异值分解的星敏感器在轨标定方法,其特征在于:步骤2.1)所述的最优星组合模型的选择是:
对于超过三颗星,有三个非零奇异值,因此组合可以由三颗星、四颗星组成;如果所有的组合都被用来构成观测模型,计算量是巨大的;此外,组合不独立,校准信息不随组合数量的增加而增加;
设视场中有n颗星星,考虑如下几种模型:
模型1:三颗星组成一个组合,如1-2-3、2-3-4、3-4-5;
模型2:
Figure FDA0004053802380000061
个星,
Figure FDA0004053802380000062
为向下取整,构成一个组合;
模型3:N-1个星构成一个组合;
模型4:与其他模型不同,该模型结合了不同数量的星形组合,如1-2-3、1-2-3-4、1-2-3-4-5;
根据可观测性分析,模型4的最小奇异值总是最大的,这意味着模型4的可观测性比其他模型好,所以将模型4作为最佳组合模型。
CN201910432037.4A 2019-05-23 2019-05-23 基于奇异值分解的星敏感器在轨标定方法 Active CN110006462B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910432037.4A CN110006462B (zh) 2019-05-23 2019-05-23 基于奇异值分解的星敏感器在轨标定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910432037.4A CN110006462B (zh) 2019-05-23 2019-05-23 基于奇异值分解的星敏感器在轨标定方法

Publications (2)

Publication Number Publication Date
CN110006462A CN110006462A (zh) 2019-07-12
CN110006462B true CN110006462B (zh) 2023-03-03

Family

ID=67177690

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910432037.4A Active CN110006462B (zh) 2019-05-23 2019-05-23 基于奇异值分解的星敏感器在轨标定方法

Country Status (1)

Country Link
CN (1) CN110006462B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113607188B (zh) * 2021-08-02 2022-07-05 北京航空航天大学 基于经纬仪叉丝成像的多视场星敏感器的标定系统及方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1923621A (zh) * 2006-10-10 2007-03-07 北京航空航天大学 基于径向排列约束的星敏感器在轨校准方法
CN103791902A (zh) * 2014-01-23 2014-05-14 中国科学院长春光学精密机械与物理研究所 适用于高机动载体的星敏感器自主导航方法
CN106403929A (zh) * 2016-08-22 2017-02-15 清华大学 星敏感器星图识别和姿态解算的快速鲁棒校验方法
CN106989761A (zh) * 2017-05-25 2017-07-28 北京航天自动控制研究所 一种基于自适应滤波的空间飞行器制导工具在轨标定方法
CN107036598A (zh) * 2017-03-30 2017-08-11 南京航空航天大学 基于陀螺误差修正的对偶四元数惯性/天文组合导航方法
CN107845096A (zh) * 2018-01-24 2018-03-27 西安平原网络科技有限公司 基于图像的行星三维信息测定方法
CN108592945A (zh) * 2018-03-27 2018-09-28 中国人民解放军国防科技大学 一种惯性/天文组合系统误差的在线标定方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1923621A (zh) * 2006-10-10 2007-03-07 北京航空航天大学 基于径向排列约束的星敏感器在轨校准方法
CN103791902A (zh) * 2014-01-23 2014-05-14 中国科学院长春光学精密机械与物理研究所 适用于高机动载体的星敏感器自主导航方法
CN106403929A (zh) * 2016-08-22 2017-02-15 清华大学 星敏感器星图识别和姿态解算的快速鲁棒校验方法
CN107036598A (zh) * 2017-03-30 2017-08-11 南京航空航天大学 基于陀螺误差修正的对偶四元数惯性/天文组合导航方法
CN106989761A (zh) * 2017-05-25 2017-07-28 北京航天自动控制研究所 一种基于自适应滤波的空间飞行器制导工具在轨标定方法
CN107845096A (zh) * 2018-01-24 2018-03-27 西安平原网络科技有限公司 基于图像的行星三维信息测定方法
CN108592945A (zh) * 2018-03-27 2018-09-28 中国人民解放军国防科技大学 一种惯性/天文组合系统误差的在线标定方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Guide Star Selection for the Three-FOV Daytime Star Sensor;Liang Wu,et al;《Sensors 2019》;20190325;第1-16页 *
On-Orbit Star Tracker Recalibration: A Case Study;Enright, J. et al;《In Proceedings of the 2015》;20150314;第1-13页 *
基于机器学习的大视场星敏感器畸变在轨标定技术;刘 源等;《红外与激光工程》;20161231;第 45 卷(第 12 期);第1-9页 *

Also Published As

Publication number Publication date
CN110006462A (zh) 2019-07-12

Similar Documents

Publication Publication Date Title
CN111156987B (zh) 基于残差补偿多速率ckf的惯性/天文组合导航方法
US9352858B2 (en) Angles-only initial orbit determination (IOD)
US9959625B2 (en) Method for fast camera pose refinement for wide area motion imagery
Gasbarri et al. Ground tests for vision based determination and control of formation flying spacecraft trajectories
CN107861919B (zh) 一种星敏感器定姿星矢量权值实时分配方法
CN107871327A (zh) 基于特征点线的单目相机位姿估计和优化方法及系统
CN113792411A (zh) 一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法
Pfeiffer et al. A computationally efficient moving horizon estimator for ultra-wideband localization on small quadrotors
CN116645392A (zh) 一种基于关键点权重的空间目标相对位姿迭代估计方法及系统
CN115560760A (zh) 一种面向无人机的视觉/激光测距高空导航方法
CN103837159A (zh) 一种经纬仪指向修正模型正交化解耦修正方法
CN110006462B (zh) 基于奇异值分解的星敏感器在轨标定方法
Chu et al. Performance comparison of tight and loose INS-Camera integration
CN111156990B (zh) 基于指向自动测定的空间碎片实时天文定位和测光方法
CN111553954A (zh) 一种基于直接法单目slam的在线光度标定方法
El-Ashmawy A comparison study between collinearity condition, coplanarity condition, and direct linear transformation (DLT) method for camera exterior orientation parameters determination
CN113916217B (zh) 一种基于分区平流层大气折射模型的星光定位方法
Stamatopoulos et al. On the self-calibration of long focal length lenses
Liu et al. Absolute navigation and positioning of mars rover using gravity-aided odometry
JP5215615B2 (ja) 3次元位置情報復元装置およびその方法
Hu et al. Initial orbit determination utilizing solution group optimization
Wishnek Optimal Information Theoretic Techniques for Electro-Optical Space Domain Awareness
Flewelling Forest inventory predictions from individual tree crowns: regression modeling within a sample framework
Liang et al. Star sensor on-orbit calibration based on multiple calibration targets
Gullu et al. Comparative analysis of least-squares approaches for 3D datum transformation in Western Turkey

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