CN112179355B - 针对光度曲线典型特征的姿态估计方法 - Google Patents

针对光度曲线典型特征的姿态估计方法 Download PDF

Info

Publication number
CN112179355B
CN112179355B CN202010910955.6A CN202010910955A CN112179355B CN 112179355 B CN112179355 B CN 112179355B CN 202010910955 A CN202010910955 A CN 202010910955A CN 112179355 B CN112179355 B CN 112179355B
Authority
CN
China
Prior art keywords
observation
station
photometric
quaternion
target
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
CN202010910955.6A
Other languages
English (en)
Other versions
CN112179355A (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.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
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 Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN202010910955.6A priority Critical patent/CN112179355B/zh
Publication of CN112179355A publication Critical patent/CN112179355A/zh
Application granted granted Critical
Publication of CN112179355B publication Critical patent/CN112179355B/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
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Astronomy & Astrophysics (AREA)
  • Automation & Control Theory (AREA)
  • General Physics & Mathematics (AREA)
  • Navigation (AREA)
  • Photometry And Measurement Of Optical Pulse Characteristics (AREA)

Abstract

本发明考虑到光度曲线中普遍存在的中断现象，以及在曲线波动上呈现大特征、中等特征和小特征的现象，本发明针对中断特征提出了“伪观测”的概念以及基于光度曲线中断互补的多站联合观测并行融合方法来解决空间目标的姿态估计问题。而当各站光度曲线不存在中断的现象时，本发明又将上述方法转化为基于多站联合观测并行融合方法，并将其用在基于大特征、中等特征和小特征光度曲线的空间目标姿态估计中。该发明相较于传统的单站光度观测方法，具有更小的姿态估计误差和更快的收敛速度，是更适合光度观测下目标姿态与角速度估计的解决方案。

Description

针对光度曲线典型特征的姿态估计方法
技术领域
本发明属于基于光学观测的空间目标姿态跟踪领域，涉及光度曲线典型特征的姿态估计方法。
背景技术
随着高轨空间目标数量的逐年增加，使得对这一区域态势感知的需求越来越高，为此需要发展相应的理论方法与技术手段。对于距离远、尺寸小的高轨空间目标，雷达受到其功率、精度的限制以及观测噪声的影响难以成为有效的观测工具。作为另一种重要手段，光学观测难以对此类目标成像，通常仅能得到反映其亮度变化的光度曲线。而光度曲线是由空间目标反射太阳光线所得到的，与目标在空间的相对位置、目标的姿态、角速度及其形状、尺寸等特征都有关系，因此光度值是随时间变化的，且在不同时段观测同一个目标所得到的光度曲线通常不同。
近年来兴起通过递推滤波技术从光度观测中实时获取目标运动信息与特征信息的研究。2009年，Wetterer等在光度观测下基于无味滤波器(Unscented Filter,UF)首次实现对火箭发动机圆柱体残骸的姿态估计；2010年，Linares等基于姿态运动学模型与轨道动力学模型通过多模型方法实现对空间目标尺寸的识别以及位置、姿态等状态量的估计，2014年又基于旋转动力学模型与轨道动力学模型通过多模型方法实现常见空间目标的形状、尺寸的识别以及状态量的估计；2015年，Holzinger等通过建立外形不确定性的一阶动力学模型利用利用粒子滤波器(Particle Filter，PF)实现了空间目标姿态快变过程中的姿态角估计；2017年，单斌等分析了典型的四棱柱、六棱柱和八棱柱空间目标光度特征和目标形状对姿态估计的影响，并探讨了算法对姿态随机缓慢机动目标的自适应跟踪能力。
从上述分析可以看出当前的研究主要集中在空间目标的建模，以及对它们的状态估计和特征识别算法上，尚未对光度曲线所具有的特征予以关注。为此本发明开展了多种因素下的目标光度曲线的分析，从中发现了光度曲线普遍具有的四种典型特征。本发明针对这些典型特征提出解决方案，该方法来自于信息融合的思想，结合光度曲线特征设计基于多站观测的并行融合估计器。
发明内容
针对光度曲线中存在的各类典型特征的技术问题，本发明提出一种针对光度曲线典型特征的姿态估计方法，该方法通过基于光度曲线中断互补的多站联合观测并行融合方法和多站联合观测并行融合方法来跟踪空间目标姿态，这两种方法的姿态估计结果收敛精度高且收敛速度快，该发明是可以实现各类典型特征光度曲线下空间目标姿态高精度估计的有效解决方案。
为实现上述目的，本发明采用以下技术方案：
针对光度曲线典型特征的姿态估计方法，包括以下步骤：
将多个观测站点联合，建立多站联合观测模型；
采用四元数对空间目标姿态进行描述，建立姿态和角速度的运动方程；
在各观测站光度曲线中断的情况下，基于多站联合观测模型、姿态和角速度的运动方程，通过伪观测建模以及基于光度曲线中断互补的多站联合观测并行融合处理，进行目标的姿态和角速度估计；
在各观测站光度曲线均不存在中断的情况下，基于多站观测的并行融合光度观测数据，进行目标的姿态和角速度估计。
作为本发明的进一步改进，多站联合观测模型过程为：
多个观测站点联合观测下同一时刻会从多个站点观测到光度数据，观测模型为：
Figure GDA0004170575680000021
其中，
Figure GDA0004170575680000022
为第N个观测站k时刻的观测数据；
Figure GDA0004170575680000023
为第N个观测站k时刻的量测噪声；
Figure GDA0004170575680000024
为第N个观测站k时刻的光度数据；光度数据计算模型为：
将每个空间目标看成是由N个平面组成，观测站测到的光度大小为各个面反射到地面观测站光度大小的总和，则光度计算模型为：
Figure GDA0004170575680000025
其中，Csun,vis＝455W/m2为可见光照射到空间目标表面单位平方米面积上的功率；Fobs(i)为太阳光线经空间目标表面i反射到地面观测站的辐射量，计算为
Figure GDA0004170575680000026
其中，Asc(i)为第i个反射面的面积；robs＝R-r为地心惯性坐标系下卫星指向地面观测站的向量，此处r和R分别是地心指向卫星的向量和地心指向观测站的向量；
Figure GDA0004170575680000027
为robs的单位向量；
Figure GDA0004170575680000031
为地心惯性坐标系下卫星表面i的单位法向量，由星体坐标系下该表面单位法向量
Figure GDA0004170575680000032
通过姿态矩阵A(q)转换得到，即
Figure GDA0004170575680000033
其中，四元数q表示目标的姿态；Fsun(i)为太阳光线辐射至卫星表面i的辐射量，计算为
Figure GDA0004170575680000034
其中，
Figure GDA0004170575680000035
表示地心惯性坐标系下卫星指向太阳的单位向量，r⊙是地心指向太阳的向量；ρtotal(i)为第i个平面的双向反射分布函数，采用Phong模型计算ρtotal，并假定光反射量主要由漫反射部分ρdiff和镜面反射部分ρspec构成：
ρtotal(i)＝ρdiff(i)+ρspec(i)。
作为本发明的进一步改进，姿态和角速度的运动方程建模过程为：
采用四元数对空间目标姿态进行描述，四元数定义为q＝[ρT q4]T，其中ρ＝[q1 q2q3]T，且满足qTq＝1的限制关系；姿态和角速度的运动方程为
Figure GDA0004170575680000036
Figure GDA0004170575680000037
其中，t为时间；
Figure GDA0004170575680000038
I是单位矩阵，对于三维向量a＝[a1 a2 a3]T，有
Figure GDA0004170575680000039
卫星角速度ω(t)＝[ωx(t) ωy(t) ωz(t)]T；w1(t)、w2(t)为零均值高斯白噪声。
作为本发明的进一步改进，还引入罗德里格斯参数GRPs，根据用到的局部误差四元数δq＝[δρT δq4]T，则GRPs表示为
Figure GDA00041705756800000310
为使误差更小，可令a＝1，f＝2(a+1)＝4；
目标姿态估计器的状态向量表示为X＝[δpT ωT]T，则姿态角与角速度联合估计的离散化动态模型为：
X(k+1)＝FX(k)+Γk
其中，X(k)为当前时刻k的目标状态向量；Γk为动态模型噪声；
Figure GDA0004170575680000041
此处
Figure GDA0004170575680000042
Δt为光度观测的采样周期，ωk为k时刻的目标角速度。
作为本发明的进一步改进，光度曲线的中断互补是利用多个观测站对空间目标进行同步观测，将每个观测站得到的光度数据分别作为量测向量中的一个观测分量；观测站选站准则为：
站点数量的选择和站点位置的设置应使光度曲线的缺失段可以互补；
由此，对于观测站{1,2,...N}，在任意时刻k，至少
Figure GDA0004170575680000043
一个观测站i使得
Figure GDA0004170575680000044
如果所加入的第二个站点未实现光度曲线中断互补，则既改变第二个观测站的位置以改变其补偿能力，或加入新的观测站，直到满足光度曲线中断互补的要求。
作为本发明的进一步改进，所述伪观测建模过程为：
在第i个观测站当前时刻k没有获得光度数据时，将该站k-1时刻对当前观测值的预测值
Figure GDA0004170575680000045
取为k时刻的观测值，并将该观测值称为伪观测建模：
Figure GDA0004170575680000046
作为本发明的进一步改进，光度曲线中断互补的并行融合步骤包括：
已知当前时刻k的n维状态向量
Figure GDA0004170575680000047
状态协方差为Pk|k；将罗德里格斯参数GRPs初始化为
Figure GDA0004170575680000048
同时将当前时刻k的四元数记为
Figure GDA0004170575680000049
则状态的sigma点由下式计算：
Figure GDA0004170575680000051
Figure GDA0004170575680000052
Figure GDA0004170575680000053
将状态的sigma点拆分为分别对应于GRPs和角速度的两部分
Figure GDA0004170575680000054
将GRPs部分的sigma点
Figure GDA0004170575680000055
转换为局部误差四元数的sigma点
Figure GDA0004170575680000056
Figure GDA0004170575680000057
则得到局部误差四元数的sigma点
Figure GDA0004170575680000058
将局部误差四元数的sigma点
Figure GDA0004170575680000059
与当前时刻四元数
Figure GDA00041705756800000510
相乘得到四元数的sigma点
Figure GDA00041705756800000511
其中，四元数乘法
Figure GDA00041705756800000512
根据上述推导，包含四元数的状态sigma点表示为
Figure GDA00041705756800000513
通过将χk(i)代入动态模型得到包含四元数的状态sigma点的一步预测：
Figure GDA00041705756800000514
令均值四元数的一步预测为
Figure GDA00041705756800000515
计算得到局部误差四元数的一步预测
Figure GDA00041705756800000516
并将其转化为GRPs
Figure GDA00041705756800000517
Figure GDA00041705756800000518
其中，四元数求逆q-1＝[-ρT q4]T；则包含GRPs的状态sigma点一步预测为
Figure GDA0004170575680000061
因此，状态一步预测均值和协方差计算为：
Figure GDA0004170575680000062
Figure GDA0004170575680000063
其中，Wi mean和Wi cov分别为计算第i采样点状态均值和协方差的权重；
若采用N个地面观测站，将包含四元数的状态sigma点一步预测χk+1(i)代入多站联合观测模型获得观测值的一步预测：
Figure GDA0004170575680000064
则N个地面站观测值的预测均值和协方差为
Figure GDA0004170575680000065
Figure GDA0004170575680000066
状态和观测值的互协方差和增益分别计算为
Figure GDA0004170575680000067
Figure GDA0004170575680000068
若k+1时刻第i个观测站的光度观测值在光度数据缺失时，伪观测值
Figure GDA0004170575680000069
在正常情况下取为实际观测值
Figure GDA00041705756800000610
则状态和状态协方差更新结果为
Figure GDA0004170575680000071
Figure GDA0004170575680000072
其中，“/”表示根据具体情况从伪观测值与实际观测值选一个作为该观测站k+1时刻的观测值。
作为本发明的进一步改进，在各观测站光度曲线均不存在中断的情况下，基于多站观测的并行融合方法，状态和状态协方差状态更新表示为
Figure GDA0004170575680000073
与现有技术相比，本发明具有以下有益效果：
本发明考虑到光度曲线中普遍存在的中断现象，以及在曲线波动上呈现大特征、中等特征和小特征的现象，本发明针对中断特征提出了“伪观测”的概念以及基于光度曲线中断互补的多站联合观测并行融合方法来解决空间目标的姿态估计问题。而当各站光度曲线不存在中断的现象时，本发明又将上述方法转化为基于多站联合观测并行融合方法，并将其用在基于大特征、中等特征和小特征光度曲线的空间目标姿态估计中。该发明相较于传统的单站光度观测方法，具有更小的姿态估计误差和更快的收敛速度，是更适合光度观测下目标姿态与角速度估计的解决方案。针对各类典型特征的光度曲线，提出基于光度曲线中断互补的多站联合观测并行融合方法和多站联合观测并行融合方法，这种途径可以保证姿态估计精度高且收敛速度快。
附图说明
图1：为光度曲线典型特征；
图2：为多站并行观测的数据处理流程图；
图3：为存在中断的光度曲线；
图4：为光度曲线存在中断时的姿态和角速度估计误差，(a)姿态估计误差；(b)角速度估计误差；
图5：为大特征光度曲线；
图6：为利用大特征光度曲线的姿态和角速度估计误差，(a)姿态估计误差；(b)角速度估计误差；
图7：为中等特征光度曲线；
图8：为利用中等特征光度曲线的姿态和角速度估计误差，(a)姿态估计误差；(b)角速度估计误差；
图9：为小特征光度曲线；
图10：为利用小特征光度曲线的姿态和角速度估计误差，(a)姿态估计误差；(b)角速度估计误差。
具体实施方式
现结合实施例、附图对本发明做进一步描述：
本发明一种针对光度曲线典型特征的姿态估计方法，步骤如下：
步骤1光度曲线典型特征的发现。
通过在不同观测时段对不同轨道倾角和升交点赤经的空间目标进行24h连续光度观测，发现空间目标光度曲线中较为普遍存在大特征、中等特征、小特征以及中断的现象，其中轨道倾角为30°的卫星的24h光度曲线如图1所示。为实现对目标姿态和角速度的有效估计，需要结合光度曲线的这些特征开展针对性的算法设计。
步骤2多站联合观测建模。
多站联合观测下同一时刻会从多个站点观测到光度数据，为便于表达及使用，给出多站联合观测模型
Figure GDA0004170575680000081
其中，
Figure GDA0004170575680000082
为第N个观测站k时刻的观测数据；
Figure GDA0004170575680000083
为第N个观测站k时刻的量测噪声；
Figure GDA0004170575680000084
为第N个观测站k时刻的光度数据。下面给出光度计算模型。
每个空间目标都可以看成是由N个平面组成的，观测站测到的光度大小为各个面反射到地面观测站光度大小的总和。光度计算模型为
Figure GDA0004170575680000085
其中，Csun,vis＝455W/m2为可见光照射到空间目标表面单位平方米面积上的功率；Fobs(i)为太阳光线经空间目标表面i反射到地面观测站的辐射量，可计算为
Figure GDA0004170575680000091
其中，Asc(i)为第i个反射面的面积；robs＝R-r为地心惯性坐标系下卫星指向地面观测站的向量，此处r和R分别是地心指向卫星的向量和地心指向观测站的向量；
Figure GDA0004170575680000092
为robs的单位向量；
Figure GDA0004170575680000093
为地心惯性坐标系下卫星表面i的单位法向量，可由星体坐标系下该表面单位法向量
Figure GDA0004170575680000094
通过姿态矩阵A(q)转换得到，即
Figure GDA0004170575680000095
其中，四元数q表示目标的姿态；Fsun(i)为太阳光线辐射至卫星表面i的辐射量，可计算为
Figure GDA0004170575680000096
其中，
Figure GDA0004170575680000097
表示地心惯性坐标系下卫星指向太阳的单位向量，r⊙是地心指向太阳的向量；ρtotal(i)为第i个平面的双向反射分布函数(BRDF)，与目标材料性质有关，本发明采用Phong模型计算ρtotal，并假定光反射量主要由漫反射部分ρdiff和镜面反射部分ρspec构成
ρtotal(i)＝ρdiff(i)+ρspec(i)
步骤3目标姿态运动学建模。
采用四元数对空间目标姿态进行描述，四元数定义为q＝[ρT q4]T，其中ρ＝[q1 q2q3]T，且满足qTq＝1的限制关系。姿态和角速度的运动方程为
Figure GDA0004170575680000098
Figure GDA0004170575680000099
其中，t为时间；
Figure GDA00041705756800000910
I是单位矩阵，对于三维向量a＝[a1 a2 a3]T，有
Figure GDA00041705756800000911
卫星角速度ω(t)＝[ωx(t) ωy(t) ωz(t)]T；w1(t)、w2(t)为零均值高斯白噪声。
为了解决四元数的乘性特性和规范化限制问题，引入广义罗德里格斯参数(GRPs)。根据UF中用到的局部误差四元数δq＝[δρT δq4]T，则GRPs表示为
Figure GDA0004170575680000101
为使误差更小，可令a＝1，f＝2(a+1)＝4。
由上述分析目标姿态估计器的状态向量可表示为X＝[δpT ωT]T，则姿态角与角速度联合估计的离散化动态模型为
X(k+1)＝FX(k)+Γk
其中，X(k)为当前时刻k的目标状态向量；Γk为动态模型噪声；
Figure GDA0004170575680000102
此处
Figure GDA0004170575680000103
Δt为光度观测的采样周期，ωk为k时刻的目标角速度。
步骤4多站观测下光度曲线的中断互补方法。
对于光度观测下较为普遍存在的光度曲线中断情况，其实质为地面观测站未接收到目标反射的太阳光线，此时光度值为无穷大。针对光度曲线中断情况，本发明提出利用多站联合观测的途径实现姿态和角速度估计，为此本发明给出选站准则。
选站准则：站点数量的选择和站点位置的设置应使光度曲线的缺失段可以互补。
由此，对于观测站{1,2,...N}，在任意时刻k，至少
Figure GDA0004170575680000105
一个观测站i使得
Figure GDA0004170575680000104
对于光度曲线中断的情况，如果所加入的第二个站点未实现光度曲线中断互补，则既可以改变第二个观测站的位置以改变其补偿能力，也可以加入新的观测站，直到满足光度曲线中断互补的要求。
本发明利用多个观测站对空间目标进行同步观测，将每个观测站得到的光度数据分别作为量测向量中的一个观测分量。
步骤5光度数据补偿的伪观测建模。
光度曲线中断时的光度数据无穷大的情况在滤波过程中会引起计算的数值问题，为此本发明提出“伪观测”概念，在第i个观测站当前时刻k没有获得光度数据时，将该站k-1时刻对当前观测值的预测值
Figure GDA0004170575680000111
取为k时刻的观测值，并将该观测值称为“伪观测”。
Figure GDA0004170575680000112
这种做法虽然使得该观测站提供的新息为零，但不会造成跟踪过程的发散，也符合实际情况。
步骤6基于光度曲线中断互补的并行融合算法。
针对光度曲线存在中断的情况，采用多站观测下光度曲线中断互补的解决方案进行目标的姿态和角速度估计。
已知当前时刻k的n维状态向量
Figure GDA0004170575680000113
状态协方差为Pk|k。将GRPs初始化为
Figure GDA0004170575680000114
同时将当前时刻k的四元数记为
Figure GDA0004170575680000115
则状态的sigma点由下式计算
Figure GDA0004170575680000116
Figure GDA0004170575680000117
Figure GDA0004170575680000118
将状态的sigma点拆分为分别对应于GRPs和角速度的两部分
Figure GDA0004170575680000119
将GRPs部分的sigma点
Figure GDA00041705756800001110
转换为局部误差四元数的sigma点
Figure GDA00041705756800001111
Figure GDA00041705756800001112
则得到局部误差四元数的sigma点
Figure GDA00041705756800001113
将局部误差四元数的sigma点
Figure GDA00041705756800001114
与当前时刻四元数
Figure GDA00041705756800001115
相乘得到四元数的sigma点
Figure GDA00041705756800001116
其中，四元数乘法
Figure GDA0004170575680000121
根据上述推导，包含四元数的状态sigma点表示为
Figure GDA0004170575680000122
通过将χk(i)代入动态模型得到包含四元数的状态sigma点的一步预测
Figure GDA0004170575680000123
令均值四元数的一步预测为
Figure GDA0004170575680000124
计算得到局部误差四元数的一步预测
Figure GDA0004170575680000125
并将其转化为GRPs
Figure GDA0004170575680000126
Figure GDA0004170575680000127
其中，四元数求逆q-1＝[-ρT q4]T。则包含GRPs的状态sigma点一步预测为
Figure GDA0004170575680000128
因此，状态一步预测均值和协方差可计算为
Figure GDA0004170575680000129
Figure GDA00041705756800001210
其中，Wi mean和Wi cov分别为计算第i采样点状态均值和协方差的权重。
若采用N个地面观测站，将包含四元数的状态sigma点一步预测χk+1(i)代入步骤2中多站联合观测模型获得观测值的一步预测
Figure GDA00041705756800001211
则N个地面站观测值的预测均值和协方差为
Figure GDA0004170575680000131
Figure GDA0004170575680000132
状态和观测值的互协方差和增益分别计算为
Figure GDA0004170575680000133
Figure GDA0004170575680000134
若k+1时刻第i个观测站的光度观测值在光度数据缺失时根据步骤4取为伪观测值
Figure GDA0004170575680000135
在正常情况下取为实际观测值
Figure GDA0004170575680000136
则状态和状态协方差更新结果为
Figure GDA0004170575680000137
Figure GDA0004170575680000138
其中的“/”表示根据具体情况从伪观测值与实际观测值选一个作为该观测站k+1时刻的观测值。
步骤7多站联合观测的并行融合算法。
在各观测站光度曲线均不存在中断的情况下，步骤5提出的方法演化为基于多站观测的并行融合方法，此时，对应步骤5中的状态更新可表示为
Figure GDA0004170575680000139
在多个传感器观测值被采用时，由于有用信息会增加，多传感器融合估计的性能要好于采用单传感器的性能。由于多站联合光度观测获得了更多的光度观测数据，采用本步骤的算法，会有助于提升空间目标姿态和角速度估计的性能。
以下结合具体实施例对本发明进行详细说明。
实施例1
光度曲线存在中断时目标姿态与角速度估计：
结合利用一段正四棱柱卫星存在中断的光度曲线对该卫星进行姿态和角速度估计的具体例子对本发明做进一步描述，实现基于光度曲线中断互补的多站联合观测并行融合方法。基本的仿真环境设定为：以倾斜地球同步轨道(IGSO)卫星为跟踪目标，卫星轨道根数设置为：半长轴a＝42166.3km，偏心率e＝0deg，倾角i＝30deg，升交点赤经Ω＝120deg，近地点幅角ω＝0deg，平近地角M＝0deg；设定卫星形状模型为正四棱柱，每个面的面积为60m2；假设每个面的镜面反射率和漫反射率都相同，分别为Rspec＝0.5和Rdiff＝0.4。卫星轨道数据由卫星工具箱(STK)仿真得到。
选择6个具有代表性的观测站对同一空间目标进行观测，为简单起见，本发明对任意第n个观测站称为站n。
选择站5在2015年5月22日05：00：00UT到2015年5月22日07：00：00UT对目标进行光度观测，发现该站观测的光度曲线存在中断，根据光度曲线互补的要求，选择站6作为第二个观测站，与站5一起进行双站联合观测。两个站点观测到的光度曲线如图3所示，从图中可以看出，每条光度曲线均有“中断”存在，且在站5的数据缺失的观测时间段内，站6可以观测到光度数据；同样，在站6的数据缺失的观测时间段内，站5也可以观测到光度数据。显然，两个观测站的观测满足互补的要求。
基于光度曲线中断互补的双站联合观测方法利用图3所示的站5和站6观测到的存在数据缺失的光度曲线进行空间目标的姿态和角速度估计。初始状态协方差P(0)＝diag(0.1,0.1,0.1,(10-5)2,(10-5)2,(10-4)2)，过程噪声协方差Q＝diag((10-4)2,(10-4)2,(10-4)2,(10-10)2,(10-10)2,(10-5)2)，每个观测站的光度观测噪声协方差为R＝0.12。
图4分别为光度数据存在缺失时的姿态和角速度估计误差，其中θ1、θ2、θ3分别为滚转、俯仰、偏航三个姿态角；ω1、ω2、ω3分别为地球惯性坐标系下x、y、z三个方向的角速度。从图中可以看出，基于站5与站6光度数据的单站观测姿态、角速度估计结果均未收敛，而基于光度曲线中断互补的双站联合观测方法确保了姿态角和角速度估计的收敛。可见在满足光度曲线中断互补的情况下，双站联合观测可以实现对姿态角与角速度的有效估计。
实施例2
大特征光度曲线下目标的姿态与角速度估计：
针对高轨空间目标，选择了基础观测站——站1，并相对于该站在不同的方位上布置其他代表性站点作为对照。最终沿经线、纬线设置了站2～5共4个观测站进行光度观测，这些站点分布在能够进行有效观测的相当大的区域范围内，而且其分布具有一定的代表性，以此对照不同观测站点间光度曲线的差异，其中站5由于中断特征在实施例1中已讨论，在本实施例中不予考虑。
选择观测区间为2015年5月22日05：00：00UT到2015年5月22日07：00：00UT时，各站观测的光度曲线如图5所示。通观图5与图3可见，如果在基准站点1和站2～5之间布置观测站点，它们的光度曲线应介于相应的站点之间，且这些站点的光度曲线相靠近；如果在比站2～5更远一些的位置上设置观测站点，则它们的光度曲线与相应的站点的观测曲线相距也不远；从空间距离上站2～5相对基准站已足够远，所以可以认为站2～5的选择在空间分布上具有代表性。同时，此时各站观测到的光度曲线均表现出共性特征：曲线波动具有明显的周期性、波峰与波谷光度差值较大。
基于双站观测方法利用图5所示的站1和站4两个代表性站点观测到的大特征光度曲线进行空间目标的姿态和角速度估计。滤波器的初始化条件与实施例1相同。
图6分别为大特征光度曲线条件下正四棱柱目标的姿态和角速度估计误差，从图中可以看出，基于站1的姿态和角速度的估计误差比站4小，其收敛速度也快一些，而双站联合观测估计误差与站1相当，但收敛速度更快。站1和站4的差别是因为光度曲线不同所导致的，这说明观测站点选址影响姿态和角速度估计精度，而双站联合观测保障了算法整体的有效性。对于实际的空间目标，观测站处于有利位置或不利位置通常是不确定的，而多站联合观测将会表现出了综合的良好性能。
对比图3和图5中光度数据的取值范围来看，光度曲线中断也属于大特征一类。针对本发明研究对象以及对应的观测时段，在关注的大范围观测区域内，各观测站光度曲线共同地表现出大特征。如果存在光度曲线中断且满足中断互补，本发明所提方案能有效估计姿态角和角速度；如果不存在光度曲线中断，对于不确定的对象，多站联合观测方法表现出综合较高的估计精度。
实施例3
中等特征光度曲线下目标姿态与角速度估计：
选择观测区间为2015年5月22日12：00：00UT到2015年5月22日14：00：00UT点时，各站观测的光度曲线如图7所示。从图中看出，此时各站观测到的光度曲线均表现出共性特征：曲线波动的周期性不是很明显、波峰与波谷光度差值不大。
基于代表性站点1和4观测到的中等特征光度曲线进行空间目标的姿态和角速度估计。在保证滤波器其他参数的初始化同实施例1的情况下，滤波器初始状态值与状态真实初始值的差值设置为实施例1滤波器的15％范围之内，可保证估计结果收敛，这说明基于中等特征光度曲线的目标姿态、角速度估计能力弱于大特征光度曲线的情况。
图8分别为中等特征光度曲线条件下的目标姿态和角速度估计误差，从图8(a)可以看出，站1和站4的单站观测姿态估计均有一个较大的偏差，而基于双站联合观测方法的姿态估计可以很快收敛。从图8(b)可以看出，双站联合观测的角速度估计效果明显好于各单站观测。由此可见，对于中等特征光度曲线，双站观测较各单站观测有更好的跟踪性能。
实施例4
小特征光度曲线下目标的姿态与角速度估计：
选择观测区间为2015年5月22日14：00：00UT到2015年5月22日16：00：00UT点时，各站观测的光度曲线如图9所示。从图中看出，此时各站观测到的光度曲线均表现出共性特征：曲线波动基本体现不出周期性、波峰与波谷光度差值很小。
基于代表性站点1和4观测到的小特征光度曲线进行空间目标的姿态和角速度估计。在保证滤波器其他参数的初始化同实施例1的情况下，滤波器初始状态值与状态真实初始值的差值设置为实施例1滤波器的5％范围之内，可保证估计结果收敛，这说明基于小特征光度曲线的目标姿态与角速度估计能力弱于大特征和中等特征光度曲线的情况。
图10为小特征光度曲线下的目标姿态和角速度估计误差，从图10(a)看出，双站观测方法的姿态角估计可以很快收敛，在偏航角和滚转角估计误差上，双站观测方法相比于单站在性能上有很大的提升。从图10(b)可以看出，双站联合观测的角速度估计误差也小于各单站观测。由此可见，对于小特征光度曲线，双站观测较各单站观测有更好的跟踪性能。
从实施例1～4的仿真结果可以看出，无论是对于光度曲线中断情况，还是大特征、中等特征和小特征光度曲线情况，本发明提出的基于并行滤波结构的多站联合光度观测方法较传统的单站光度观测方法均表现出了更小的姿态和角速度估计误差和更快的收敛速度。此外，中等特征光度曲线和小特征光度曲线情况下滤波器初始状态值与状态真实初始值的差值只有分别设置为大特征光度曲线情况的15％和5％时，才能保证姿态和角速度估计结果收敛，由此可以看出，随着目标光度曲线特征的减弱，滤波器对目标姿态和角速度估计的能力减弱。
对于多种因素影响下难以预测的光度曲线中断现象，本发明提出的基于光度曲线中断互补多站联合观测的并行融合方法是有效的解决方案。对基于光度观测的高轨空间目标的姿态和角速度估计，应该优先选择大特征的目标光度曲线，这将提升滤波器对目标未知状态估计的能力。在目标的光度曲线为小特征或中等特征时，采用本发明所提出的基于多站联合观测的并行融合方法也可以实现对目标姿态的有效估计。
以上披露的所有文章和参考资料，包括专利申请和出版物，出于各种目的通过援引结合于此。描述组合的术语“基本由…构成”应该包括所确定的元件、成分、部件或步骤以及实质上没有影响该组合的基本新颖特征的其他元件、成分、部件或步骤。使用术语“包含”或“包括”来描述这里的元件、成分、部件或步骤的组合也想到了基本由这些元件、成分、部件或步骤构成的实施方式。这里通过使用术语“可以”，旨在说明“可以”包括的所描述的任何属性都是可选的。
多个元件、成分、部件或步骤能够由单个集成元件、成分、部件或步骤来提供。另选地，单个集成元件、成分、部件或步骤可以被分成分离的多个元件、成分、部件或步骤。用来描述元件、成分、部件或步骤的公开“一”或“一个”并不说为了排除其他的元件、成分、部件或步骤。
应该理解，以上描述是为了进行图示说明而不是为了进行限制。通过阅读上述描述，在所提供的示例之外的许多实施例和许多应用对本领域技术人员来说都将是显而易见的。因此，本教导的范围不应该参照上述描述来确定，而是应该参照前述权利要求以及这些权利要求所拥有的等价物的全部范围来确定。出于全面之目的，所有文章和参考包括专利申请和公告的公开都通过参考结合在本文中。在前述权利要求中省略这里公开的主题的任何方面并不是为了放弃该主体内容，也不应该认为申请人没有将该主题考虑为所公开的发明主题的一部分。

Claims (4)

1.针对光度曲线典型特征的姿态估计方法，其特征在于，包括以下步骤：
将多个观测站点联合，建立多站联合观测模型；
采用四元数对空间目标姿态进行描述，建立姿态和角速度的运动方程；
在各观测站光度曲线中断的情况下，基于多站联合观测模型、姿态和角速度的运动方程，通过伪观测建模以及基于光度曲线中断互补的多站联合观测并行融合处理，进行目标的姿态和角速度估计；
光度曲线中断互补的并行融合步骤包括：
已知当前时刻k的n维状态向量
Figure FDA0004113081650000011
其中
Figure FDA0004113081650000012
为用广义罗德里格斯参数GRPs表示的姿态角，ωk|k为角速度，状态协方差为Pk|k；将GRPs初始化为
Figure FDA0004113081650000013
同时将当前时刻k的姿态四元数记为
Figure FDA0004113081650000014
则状态的sigma点由下式计算：
Figure FDA0004113081650000015
Figure FDA0004113081650000016
Figure FDA0004113081650000017
将状态的sigma点拆分为分别对应于GRPs和角速度的两部分
Figure FDA0004113081650000018
将GRPs部分的sigma点
Figure FDA0004113081650000019
转换为局部误差四元数的sigma点
Figure FDA00041130816500000110
Figure FDA00041130816500000111
则得到局部误差四元数的sigma点
Figure FDA00041130816500000112
将局部误差四元数的sigma点
Figure FDA00041130816500000113
与当前时刻四元数
Figure FDA00041130816500000114
相乘得到四元数的sigma点
Figure FDA00041130816500000115
其中，四元数乘法
Figure FDA00041130816500000116
根据上述推导，包含四元数的状态sigma点表示为
Figure FDA0004113081650000021
通过将χk(i)代入动态模型得到包含四元数的状态sigma点的一步预测：
Figure FDA0004113081650000022
令均值四元数的一步预测为
Figure FDA0004113081650000023
计算得到局部误差四元数的一步预测
Figure FDA0004113081650000024
并将其转化为GRPs
Figure FDA0004113081650000025
Figure FDA0004113081650000026
其中，四元数q＝[ρT q4]T求逆q-1＝[-ρT q4]T；则包含GRPs的状态sigma点一步预测为
Figure FDA0004113081650000027
因此，状态一步预测均值和协方差计算为：
Figure FDA0004113081650000028
Figure FDA0004113081650000029
其中，Q为过程噪声协方差，Wi mean和Wi cov分别为计算第i采样点状态均值和协方差的权重；
若采用N个地面观测站，将包含四元数的状态sigma点一步预测χk+1(i)代入多站联合观测模型获得观测值的一步预测：
Figure FDA00041130816500000210
则N个地面站观测值的预测均值和协方差为
Figure FDA0004113081650000031
Figure FDA0004113081650000032
其中，R为观测噪声协方差；
状态和观测值的互协方差和增益分别计算为
Figure FDA0004113081650000033
Figure FDA0004113081650000034
若k+1时刻第i个观测站的光度观测值在光度数据缺失时，伪观测值
Figure FDA0004113081650000035
在正常情况下取为实际观测值
Figure FDA0004113081650000036
则状态和状态协方差更新结果为
Figure FDA0004113081650000037
Figure FDA0004113081650000038
其中，“/”表示根据具体情况从伪观测值与实际观测值选一个作为该观测站k+1时刻的观测值；
在各观测站光度曲线均不存在中断的情况下，基于多站观测的并行融合光度观测数据，进行目标的姿态和角速度估计；
建立多站联合观测模型过程为：
多个观测站点联合观测下同一时刻会从多个站点观测到光度数据，观测模型为：
Figure FDA0004113081650000039
其中，
Figure FDA00041130816500000310
为第i个观测站k时刻的观测数据,i∈{1,…,N},N为观测站点的数量；
Figure FDA00041130816500000311
为第i个观测站k时刻的量测噪声；
Figure FDA0004113081650000041
为第i个观测站k时刻的光度数据；光度数据计算模型为：
将每个空间目标看成是由M个平面组成，观测站测到的光度大小为各个面反射到地面观测站光度大小的总和，则光度计算模型为：
Figure FDA0004113081650000042
其中，Csun,vis＝455W/m2为可见光照射到空间目标表面单位平方米面积上的功率；Fobs(i)为太阳光线经空间目标表面i反射到地面观测站的辐射量，计算为
Figure FDA0004113081650000043
其中，Asc(i)为第i个反射面的面积；robs＝R-r为地心惯性坐标系下卫星指向地面观测站的向量，此处r和R分别是地心指向卫星的向量和地心指向观测站的向量；
Figure FDA0004113081650000044
为robs的单位向量；
Figure FDA0004113081650000045
为地心惯性坐标系下卫星表面i的单位法向量，由星体坐标系下该表面单位法向量
Figure FDA0004113081650000046
通过姿态矩阵A(q)转换得到，即
Figure FDA0004113081650000047
其中，四元数q表示目标的姿态；Fsun(i)为太阳光线辐射至卫星表面i的辐射量，计算为
Figure FDA0004113081650000048
其中，
Figure FDA0004113081650000049
表示地心惯性坐标系下卫星指向太阳的单位向量，r⊙是地心指向太阳的向量；ρtotal(i)为第i个平面的双向反射分布函数，采用Phong模型计算ρtotal，并假定光反射量主要由漫反射部分ρdiff和镜面反射部分ρspec构成：
ρtotal(i)＝ρdiff(i)+ρspec(i)；
姿态和角速度的运动方程建模过程为：
采用四元数对空间目标姿态进行描述，四元数定义为q＝[ρT q4]T，其中ρ＝[q1 q2 q3]T，且满足qTq＝1的限制关系；姿态和角速度的运动方程为
Figure FDA0004113081650000051
Figure FDA0004113081650000052
其中，t为时间；
Figure FDA0004113081650000053
I是单位矩阵；卫星角速度ω(t)＝[ωx(t) ωy(t) ωz(t)]T；w1(t)、w2(t)为零均值高斯白噪声；
还引入罗德里格斯参数GRPs，根据用到的局部误差四元数δq＝[δρT δq4]T，则GRPs表示为
Figure FDA0004113081650000054
为使误差更小，可令a＝1，f＝2(a+1)＝4；
目标姿态估计器的状态向量表示为X＝[δpT ωT]T，则姿态角与角速度联合估计的离散化动态模型为：
X(k+1)＝FX(k)+Γk
其中，X(k)为当前时刻k的目标状态向量；Γk为动态模型噪声；
Figure FDA0004113081650000055
此处
Figure FDA0004113081650000056
Δt为光度观测的采样周期，ωk为k时刻的目标角速度。
2.根据权利要求1所述的针对光度曲线典型特征的姿态估计方法，其特征在于，
光度曲线的中断互补是利用多个观测站对空间目标进行同步观测，将每个观测站得到的光度数据分别作为量测向量中的一个观测分量；观测站选站准则为：
站点数量的选择和站点位置的设置应使光度曲线的缺失段可以互补；
由此，对于观测站{1,2,...N}，在任意时刻k，至少一个观测站i使得
Figure FDA0004113081650000057
如果所加入的第二个站点未实现光度曲线中断互补，则改变第二个观测站的位置以改变其补偿能力，或加入新的观测站，直到满足光度曲线中断互补的要求。
3.根据权利要求2所述的针对光度曲线典型特征的姿态估计方法，其特征在于，所述伪观测建模过程为：
在第i个观测站当前时刻k没有获得光度数据时，将该站k-1时刻对当前观测值的预测值
Figure FDA0004113081650000061
取为k时刻的观测值，并将该观测值称为伪观测建模：
Figure FDA0004113081650000062
4.根据权利要求1所述的针对光度曲线典型特征的姿态估计方法，其特征在于，在各观测站光度曲线均不存在中断的情况下，基于多站观测的并行融合方法，状态更新表示为
Figure FDA0004113081650000063
CN202010910955.6A 2020-09-02 2020-09-02 针对光度曲线典型特征的姿态估计方法 Active CN112179355B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010910955.6A CN112179355B (zh) 2020-09-02 2020-09-02 针对光度曲线典型特征的姿态估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010910955.6A CN112179355B (zh) 2020-09-02 2020-09-02 针对光度曲线典型特征的姿态估计方法

Publications (2)

Publication Number Publication Date
CN112179355A CN112179355A (zh) 2021-01-05
CN112179355B true CN112179355B (zh) 2023-05-26

Family

ID=73925571

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010910955.6A Active CN112179355B (zh) 2020-09-02 2020-09-02 针对光度曲线典型特征的姿态估计方法

Country Status (1)

Country Link
CN (1) CN112179355B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112926237B (zh) * 2021-01-28 2024-05-24 南京航空航天大学 一种基于光度信号的空间目标关键特征辨识方法
WO2024116361A1 (ja) * 2022-12-01 2024-06-06 三菱電機株式会社 形状姿勢推定装置、形状姿勢推定方法及び形状姿勢推定システム

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108415098A (zh) * 2018-02-28 2018-08-17 西安交通大学 基于光度曲线对空间高轨小尺寸目标特征识别方法

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6166810A (en) * 1997-12-05 2000-12-26 Nippon Telegraph And Telephone Corporation Method and apparatus for determining distance
FR2902556B1 (fr) * 2006-06-15 2008-08-15 Valeo Vision Sa Procede de determination d'une distance de visibilite pour un conducteur d'un vehicule.
US8222582B1 (en) * 2008-09-30 2012-07-17 Anderson Mark J Celestial navigation using stellar narrow-band emission
ES2399636T3 (es) * 2009-07-29 2013-04-02 Metaio Gmbh Método para determinar la postura de una cámara con respecto a por lo menos un objeto real
CN105354875B (zh) * 2015-09-25 2018-01-23 厦门大学 一种室内环境二维与三维联合模型的构建方法和系统
JP6761244B2 (ja) * 2015-12-28 2020-09-23 川崎重工業株式会社 乗物
US10121248B2 (en) * 2016-08-19 2018-11-06 Raytheon BBN Technologies, Corp. Automated system and method for determining positional order through photometric and geospatial data
JP6699567B2 (ja) * 2017-01-17 2020-05-27 トヨタ自動車株式会社 撮像装置
CN108253962A (zh) * 2017-12-18 2018-07-06 中北智杰科技（北京）有限公司 一种低照度环境下新能源无人驾驶汽车定位方法
CN110618466B (zh) * 2018-06-20 2021-06-18 天津工业大学 一种空间目标姿态可探测性度量方法
CN109633724B (zh) * 2019-01-16 2023-03-03 电子科技大学 基于单星与多地面站联合测量的无源目标定位方法
CN111507132B (zh) * 2019-01-31 2023-07-07 杭州海康机器人股份有限公司 一种定位方法、装置及设备
CN110059292B (zh) * 2019-04-24 2023-01-24 中国人民解放军战略支援部队航天工程大学 一种空间目标姿态识别方法
CN110146082B (zh) * 2019-05-05 2021-03-19 中国人民解放军63921部队 利用测速数据实时估计航天器异常姿态的方法和设备

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108415098A (zh) * 2018-02-28 2018-08-17 西安交通大学 基于光度曲线对空间高轨小尺寸目标特征识别方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
单斌 ; 梁勇奇 ; 李恒年 ; .基于光度观测的空间目标姿态与角速度估计.光学学报.2017,37(05),全文. *

Also Published As

Publication number Publication date
CN112179355A (zh) 2021-01-05

Similar Documents

Publication Publication Date Title
Opromolla et al. Pose estimation for spacecraft relative navigation using model-based algorithms
CN106338753B (zh) 一种基于地面站/星间链路/gnss联合测量的地球同步轨道星座定轨方法
US8494687B2 (en) Method for enhancing a three dimensional image from a plurality of frames of flash LIDAR data
Trebi-Ollennu et al. Design and analysis of a sun sensor for planetary rover absolute heading detection
CN106595674B (zh) 基于星敏感器和星间链路的heo卫星编队飞行自主导航方法
CN111102981B (zh) 一种基于ukf的高精度卫星相对导航方法
CN112179355B (zh) 针对光度曲线典型特征的姿态估计方法
Pittet et al. Spin motion determination of the Envisat satellite through laser ranging measurements from a single pass measured by a single station
CN108827322B (zh) 一种多星协同测向定位观测系统优化设计与评估方法
Ely et al. Altair navigation during translunar cruise, lunar orbit, descent, and landing
CN115790575B (zh) 一种基于多星协同无源探测的巨型星座目标跟踪方法
CN111125874A (zh) 一种动平台高精度测轨预报方法
Dietrich et al. Orbit determination using flash lidar around small bodies
CN110146092B (zh) 基于导航信息评价的双体小行星探测轨迹优化方法
CN115184916A (zh) 一种海面风速联合反演方法、装置、介质及计算设备
CN101813481A (zh) 用于机载的基于虚拟水平基准修正的惯性与天文定位方法
US12078716B2 (en) System and method of hypersonic object tracking
Witte et al. No GPS? No problem! Exploring the dunes of titan with dragonfly using visual odometry
CN116698048A (zh) 一种基于脉冲星/星间测距/陆标的组合导航方法
CN114964215A (zh) 多个探测目标主体的轨道确定方法及装置、电子设备
Jovanovic et al. MISR photogrammetric data reduction for geophysical retrievals
Shuang et al. Autonomous optical navigation for landing on asteroids
Chee et al. Norm-constrained unscented kalman filter with application to high area-to-mass ratio space-debris tracking
Zhai et al. Role of Topocentric Parallax in Near-Earth Object Initial Orbit Determination
CN112926237B (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