CN110017832B - 一种基于Gauss解群优选的短弧初轨确定方法 - Google Patents

一种基于Gauss解群优选的短弧初轨确定方法 Download PDF

Info

Publication number
CN110017832B
CN110017832B CN201910207587.6A CN201910207587A CN110017832B CN 110017832 B CN110017832 B CN 110017832B CN 201910207587 A CN201910207587 A CN 201910207587A CN 110017832 B CN110017832 B CN 110017832B
Authority
CN
China
Prior art keywords
track
data
dimensional
group
component vector
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
CN201910207587.6A
Other languages
English (en)
Other versions
CN110017832A (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.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and 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 Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN201910207587.6A priority Critical patent/CN110017832B/zh
Publication of CN110017832A publication Critical patent/CN110017832A/zh
Priority to US16/762,505 priority patent/US11662209B2/en
Priority to PCT/CN2019/107186 priority patent/WO2020186719A1/zh
Application granted granted Critical
Publication of CN110017832B publication Critical patent/CN110017832B/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/20—Instruments for performing navigational calculations
    • 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/02—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by astronomical means
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/242—Orbits and trajectories
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G3/00—Observing or tracking cosmonautic vehicles
    • 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

Abstract

本发明公开一种基于Gauss解群优选的短弧初轨确定方法，属于天文动力学领域。包括：将观测数据进行分组，对每组数据，采用Gauss方法求出对应时刻的目标状态矢量，构成初步估计的解集；将初步估计的解集拆分为位置分矢量解集和速度分矢量解集，分别进行分群，得到位置分矢量解群和速度分矢量解群；基于位置分矢量解群和速度分矢量解群，生成二维轨迹解集；采用轨迹优选的方法对每条二维轨迹进行评估，计算最优二维轨迹对应轨道根数，完成初轨确定。本发明对观测数据进行分组，对每组数据调用Gauss，加大数据的利用程度，结合多组数据最终确定的短弧目标轨道根数更贴近真实。采用Gauss解群优选进行轨迹评估，最优二维轨迹的轨道根数作为估计解，解决了多根问题。

Description

一种基于Gauss解群优选的短弧初轨确定方法
技术领域
本发明属于天文动力学领域，更具体地，涉及一种基于Gauss解群优选的短弧初轨确定方法。
背景技术
随着天文动力学的发展，初轨确定技术已成为当今最重要的技术之一。传统的初轨确定方法是通过地面观测数据计算卫星或空间碎片的轨道，这些方法通常是基于两体问题的运动方程。由于天基观测系统观测时刻短、可获取的数据较少，导致传统的初轨确定方法不适用于天基观测数据；此外，如果待检测轨迹的曲率过小，则有可能无法确定初轨。
当前的初轨确定方法主要包括：传统的Laplace法、Gauss法及相关改进算法。传统Gauss法的缺陷是：求解高阶多项式后所得的根太多，需要进行筛选；改进Gauss法可以避免重根问题，但需要选取合适的初值，否则会使得算法迭代发散或最终取的平凡解甚至错误解。
当前研究资料表明，在基于地面观测数据的初轨确定领域，已有许多性能优异的方法，而在基于天基观测数据的初轨确定领域，相关研究算法的时间间隔要求都比较长，对于短弧初轨确定适用性不高，本方法可解决传统Gauss方法中存在的多根问题，并适用于短弧初轨确定任务。
发明内容
针对现有技术的缺陷，本发明的目的在于解决现有技术基于传统Gauss方法的短弧初轨确定方法误差较大、难以从多组解中有效选择的技术问题。
为实现上述目的，第一方面，本发明实施例提供了一种基于Gauss解群优选的短弧初轨确定方法，该方法包括以下步骤：
S1.将观测数据进行分组，对每组数据，采用Gauss方法求出对应时刻的目标状态矢量，构成初步估计的解集；
S2.将初步估计的解集拆分为位置分矢量解集和速度分矢量解集，分别进行分群，得到位置分矢量解群和速度分矢量解群；
S3.基于位置分矢量解群和速度分矢量解群，生成二维轨迹解集；
S4.采用轨迹优选的方法对每条二维轨迹进行评估，计算最优二维轨迹对应轨道根数，完成初轨确定。
具体地，步骤S1包括以下子步骤：
S101.对观测数据进行分组，每3个时刻对应的矢量数据分为一组；
S102.对每组矢量数据，采用传统Gauss法求出对应时刻的目标状态矢量，构成初步估计的解集。
具体地，假定观测时刻为[t0,tk]，k+1≥3，对于[t0,tk]共k+1个时刻的输入数据，取间隔参数
Figure BDA0001999469440000026
其中，
Figure BDA0001999469440000027
表示向下取整，rate为间隔率，分组后的观测数据表示为
Figure BDA0001999469440000021
其中，
Figure BDA0001999469440000022
表示卫星观测位置矢量，
Figure BDA0001999469440000023
表示观测角度矢量，
Figure BDA0001999469440000024
Figure BDA0001999469440000025
具体地，步骤S2包括以下子步骤：
S201.剔除初步估计的解集中不合理的解；
S202.将剩余的状态矢量解集拆分为位置分矢量解集和速度分矢量解集，分别进行分群，使得正确的解尽可能在同一个群中，得到位置分矢量解群和速度分矢量解群；
S203.对于位置分矢量解群和速度分矢量解群分别进行降噪预处理。
具体地，剔除满足
Figure BDA0001999469440000031
或者
Figure BDA0001999469440000032
或者
Figure BDA0001999469440000033
任一条件的位置数据和速度数据，其中，
Figure BDA0001999469440000034
表示ts时刻的目标位置矢量，
Figure BDA0001999469440000035
表示ts时刻的目标速度矢量。
具体地，步骤S202中采用k-means聚类方法进行聚类分群，并采用Chauvenet's-criterion判别进行异常数据消除。
具体地，步骤S3包括以下子步骤：
S301.对降噪后的位置分矢量解群和速度分矢量解群进行拟合，构建各个时刻的状态矢量组合；
S302.根据各个时刻对应的状态矢量组合，生成目标轨道三维轨迹解集；
S303.根据观测平台测量状态，将目标轨道三维轨迹解集投影到瞬时观测像面，得到二维轨迹解集。
具体地，步骤S4包括以下子步骤：
S401.对每条二维轨迹,计算导数误差与位置误差；
第m条二维轨迹的导数误差的计算公式如下：
Figure BDA0001999469440000036
Figure BDA0001999469440000037
Figure BDA0001999469440000038
Figure BDA0001999469440000039
Figure BDA00019994694400000310
其中，U″*，V″*，U'*，V'*为二维轨迹
Figure BDA00019994694400000311
的二次多项式拟合参数，
Figure BDA00019994694400000312
为估计轨迹l′m在tn时刻对应的像面坐标，m表示第m条二维轨迹，k为观测时刻总数；
第m条二维轨迹的位置误差的计算公式如下：
Figure BDA0001999469440000041
Figure BDA0001999469440000042
Figure BDA0001999469440000043
其中，
Figure BDA0001999469440000044
为预测值与观测值之间的像面距离，
Figure BDA0001999469440000045
为估计轨迹与观测轨迹交点和观测点之间的像面距离，tn为观测时刻，
Figure BDA0001999469440000046
为估计轨迹l′m在tn时刻对应的像面坐标，
Figure BDA0001999469440000047
为观测轨迹l*在tn时刻对应的像面坐标；
Figure BDA0001999469440000048
为估计轨迹l′m与观测轨迹交点的像面坐标；
S402.综合考虑导数误差与位置误差，选择最优二维轨迹，将其对应的轨道根数作为最佳估计值；
S403.输出最优二维轨迹对应轨道根数，完成初轨确定。
具体地，步骤S402具体为：
对于第m条轨道，分别对导数误差和位置误差从小到大进行排序，得到排名序号Rankvel,m和Rankdis,m，两者在1～Nsv之间，并计算两者之和作为最终误差评价：
Rankm＝Rankvel,m+Rankdis,m
对于最终误差评价，取其最小者mopt作为最佳预测轨迹：
Figure BDA0001999469440000049
其中，Nsv为各个时刻对应的目标位置、速度数据组合所求得的轨道根数的个数。
第二方面，本发明实施例提供了一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器执行时实现上述第一方面所述的短弧初轨确定方法。
总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：
1.本发明通过对观测数据进行分组，并对每组数据调用传统Gauss方法，加大数据的利用程度，结合多组数据最终确定的短弧目标轨道根数会更加贴近真实值。
2.本发明针对多根问题，将所有根经过剔除不合理数据、降噪平滑等处理方式后，计算它们对应的轨道根数，并采用轨迹优先的方法去进行评估，最终选择评价指标最高的轨道根数作为估计解，这一做法能够有效地解决传统Gauss法中存在的多根问题。
附图说明
图1为本发明实施例提供的一种基于Gauss解群优选的短弧初轨确定方法流程图；
图2为本发明实施例提供的天基观测示意图；
图3为本发明实施例提供的Gauss法解得的目标位置解群图；
图4为本发明实施例提供的位置解群进行k均值聚类后的聚类图；
图5为本发明实施例提供的采用Savitzky-Golay平滑后的位置解群图；
图6为本发明实施例提供的计算轨道根数得到各位置解对应的轨道三维图；
图7为本发明实施例提供的各位置解对应三维轨道投影到像面上的二维图；
图8为本发明实施例提供的最优轨道与实际观测轨道的三维图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
首先，对本发明涉及的一些术语进行解释。
解群：解的集合。
如图1所示，本发明提供了一种基于Gauss解群优选的短弧初轨确定方法，该方法包括以下步骤：
S1.将观测数据进行分组，对每组数据，采用Gauss方法求出对应时刻的目标状态矢量，构成初步估计的解集；
S2.将初步估计的解集拆分为位置分矢量解集和速度分矢量解集，分别进行分群，得到位置分矢量解群和速度分矢量解群；
S3.基于位置分矢量解群和速度分矢量解群，生成二维轨迹解集；
S4.采用轨迹优选的方法对每条二维轨迹进行评估，计算最优二维轨迹对应轨道根数，完成初轨确定。
步骤S1.将观测数据进行分组，对每组数据，采用Gauss方法求出对应时刻的目标状态矢量，构成初步估计的解集。
S101.对观测数据进行分组，每3个时刻对应的矢量数据分为一组。
观测数据是天基单目成像平台观测数据或者地面观测数据。所述观测数据包括：各个观测时刻对应的观测位置(平台本身相对地心坐标系位置)、观测角度矢量，分别用变量
Figure BDA0001999469440000064
表示ts时刻卫星观测位置矢量、
Figure BDA0001999469440000065
表示ts时刻观测角度矢量。同一时刻下观测卫星的位置、观测相机的角度以及观测成像数据共同构成一对数据，多个时刻数据构成一组数据。
为了后续构建解群，本发明首先需要对观测数据进行分组，每3个时刻对应的矢量数据分为一组。假定观测时刻为[t0,tk](k+1≥3)，对于[t0,tk]共k+1个时刻的输入数据，取间隔参数
Figure BDA0001999469440000066
其中，
Figure BDA0001999469440000067
表示向下取整，rate为间隔率，通常取值0.2。分组后的观测数据表示为
Figure BDA0001999469440000061
其中，
Figure BDA0001999469440000062
表示卫星观测位置矢量，
Figure BDA0001999469440000063
表示观测角度矢量，
Figure BDA0001999469440000071
Figure BDA0001999469440000072
S102.对每组矢量数据，采用传统Gauss法求出对应时刻的目标状态矢量，构成初步估计的解集。
由于传统Gauss方法输入仅为3个数据，所输入的信息有限，本发明考虑多次使用Gauss方法以加大数据的利用程度，结合多组数据最终确定的短弧目标轨道根数会更加贴近真实值。
Gauss法是通过输入三个不同时刻对应的观测位置向量和观测角度向量，求解中间时刻s2的目标位置矢量与目标速度矢量。目标指无动力目标，可以是碎片或者失控卫星等。按照本发明方法将一系列时刻中每三个分为一组，每组数据均采用一次传统Gauss法求解得到该组数据中间时刻的目标位置矢量和速度矢量解。每组矢量数据得到的解的个数可能为1、2、3。多次调用Gauss方法，求出不同时刻对应的目标矢量构成初步估计的解的集合。
步骤S2.将初步估计的解集拆分为位置分矢量解集和速度分矢量解集，分别进行分群，得到位置分矢量解群和速度分矢量解群。
S201.剔除初步估计的解集中不合理的解。
剔除满足
Figure BDA0001999469440000073
(地球半径)或者
Figure BDA0001999469440000074
(地球同步静止轨道半径)或者
Figure BDA0001999469440000075
(第二宇宙速度)任一条件的位置数据和速度数据，其中，
Figure BDA0001999469440000076
表示ts时刻的目标位置矢量，
Figure BDA0001999469440000077
表示ts时刻的目标速度矢量。
很显然，满足上述三个条件任一个的解，是不符合常理的解。因此，剔除所求目标位置距离地心小于地球半径或大于地球同步静止轨道半径的目标位置数据以及大于第二宇宙速度的目标速度数据。
S202.将剩余的状态矢量解集拆分为位置分矢量解集和速度分矢量解集，分别进行分群，使得正确的解尽可能在同一个群中，得到位置分矢量解群和速度分矢量解群。
为了将正确的根和错误的根分开，尽量使得正确的根在同一个群，对位置分矢量解集和速度分矢量解集分别进行分群。本发明优选k-means聚类方法进行聚类分群(k取值2或者3)，并采用Chauvenet's-criterion(肖维奈准则)判别进行异常数据消除，异常数据是指同一个群中表现异常的数据，例如，离群数据。
若每个时刻对应的位置解个数小于等于1，则不需要进行k均值聚类；若其中任意时刻对应的位置解个数最大为2，则进行一次k＝2的k均值聚类；若其中任意时刻对应的位置解个数最大为3，则进行一次k＝3的k均值聚类。速度分矢量解集聚类同上。
S203.对于位置分矢量解群和速度分矢量解群分别进行降噪预处理。
本发明优选基于Savitzky-Golay滤波方法进行降噪。
所有根经过剔除不合理数据、降噪平滑等处理方式后，能够提升轨迹优选的精确度。
步骤S3.基于位置分矢量解群和速度分矢量解群，生成二维轨迹解集。
S301.对降噪后的位置分矢量解群和速度分矢量解群进行拟合，构建各个时刻的状态矢量组合。
利用降噪后的数据，进行目标位置分矢量拟合以及目标速度分矢量拟合。时刻ts对应的速度矢量和位置矢量组合为
Figure BDA0001999469440000081
S302.根据各个时刻对应的状态矢量组合，生成目标轨道三维轨迹解集。
通过某一时刻ts对应的速度矢量和位置矢量组合
Figure BDA0001999469440000082
计算出
Figure BDA0001999469440000083
对应的轨道根数
Figure BDA0001999469440000084
通过轨道根数可以确定三维轨迹解集
Figure BDA0001999469440000085
其中，a为半长轴，e为第一偏心率，i为轨道倾角，Ω为升交点赤经，ω为近地点幅角，M为平近点角，
Figure BDA0001999469440000091
为观测轨迹
Figure BDA0001999469440000092
在tn时刻对应的三维坐标，m＝1,2,…,Nsv，n＝0,1,…,k，Nsv为各个时刻对应的目标位置、速度数据组合所求得的轨道根数的总个数。
S303.根据观测平台测量状态，将目标轨道三维轨迹解集投影到瞬时观测像面，得到二维轨迹解集。
目标轨道三维轨迹解集
Figure BDA0001999469440000093
投影到瞬时观测像面上，即可得到二维轨迹解集
Figure BDA0001999469440000094
其中，
Figure BDA0001999469440000095
为第m条三维观测轨迹，
Figure BDA0001999469440000096
为观测轨迹
Figure BDA0001999469440000097
在在tn时刻对应的三维坐标，l′m为第m条估计轨迹，
Figure BDA0001999469440000098
为估计轨迹l′m在tn时刻对应的像面坐标，n＝0,1,…,k，m＝1,2,…,Nsv。
步骤S4.采用轨迹优选的方法对每条二维轨迹进行评估，计算最优二维轨迹对应轨道根数，完成初轨确定。
S401.对每条二维轨迹,计算导数误差与位置误差。
第m条二维轨迹的导数误差的计算公式如下：
Figure BDA0001999469440000099
Figure BDA00019994694400000910
Figure BDA00019994694400000911
Figure BDA00019994694400000912
Figure BDA00019994694400000913
其中，U″*，V″*，U'*，V'*为二维轨迹
Figure BDA00019994694400000914
的二次多项式拟合参数，
Figure BDA00019994694400000915
为估计轨迹l′m在tn时刻对应的像面坐标，m表示第m条二维轨迹，k为观测时刻总数。
相机参数(朝向)保持不变，可以得到各个时刻对目标观测的相面坐标
Figure BDA0001999469440000101
对于这些观测数据在U,V方向分别做二次多项式拟合，可以求得第二项、第四项。
第m条二维轨迹的位置误差的计算公式如下：
Figure BDA0001999469440000102
Figure BDA0001999469440000103
Figure BDA0001999469440000104
其中，
Figure BDA0001999469440000105
为预测值与观测值之间的像面距离，
Figure BDA0001999469440000106
为估计轨迹与观测轨迹交点和观测点之间的像面距离，tn为观测时刻，
Figure BDA0001999469440000107
为估计轨迹l′m在tn时刻对应的像面坐标，
Figure BDA0001999469440000108
为观测轨迹l*在tn时刻对应的像面坐标，
Figure BDA0001999469440000109
为估计轨迹l′m与观测轨迹交点的像面坐标。
S402.综合考虑导数误差与位置误差，选择最优二维轨迹，将其对应的轨道根数作为最佳估计值。
对于第m条轨道，分别对导数误差和位置误差从小到大进行排序，得到排名序号Rankvel,m和Rankdis,m，两者在1～Nsv之间，并计算两者之和作为最终误差评价：
Rankm＝Rankvel,m+Rankdis,m
对于最终误差评价，取其最小者mopt作为最优二维轨迹：
Figure BDA00019994694400001010
其中，Nsv为各个时刻对应的目标位置、速度数据组合所求得的轨道根数的个数。
S403.输出最优二维轨迹对应轨道根数，完成初轨确定。
适用于单一天基成像观测平台，短时间观测条件(大于5s)下的远距离无动力目标，如碎片、失控卫星等的初轨确定。本发明能够有效解决传统Gauss方法难以解决的多根问题，对于原始观测数据的利用更加充分，利用解群的概念和优选方法提高了初轨定轨精度，在相同输入误差条件下，本技术方案对短时间观测数据比其他方案表现更佳，即本发明较其他发明对于短弧观测更适用。
实施例1
如表1所示，为实施例1中观测卫星与目标轨道根数数据。
表1
Figure BDA0001999469440000111
S101.对观测数据进行分组，每3个时刻对应的矢量数据分为一组。
对于[t0,t170]共171个时刻的输入数据，取间隔参数Nrate＝34，将观测数据
Figure BDA0001999469440000112
分为一组，其中，
Figure BDA0001999469440000113
表示向下取整，rate为间隔率取rate＝0.2，
Figure BDA0001999469440000114
Figure BDA0001999469440000115
表示卫星观测位置向量，
Figure BDA0001999469440000116
表示观测角度向量。
对于表1数据下的观测卫星与目标轨道根数，仿真得到表2中观测卫星与目标部分位置数据。
表2
Figure BDA0001999469440000117
Figure BDA0001999469440000121
步骤S102.对每组矢量数据，采用传统Gauss法求出对应时刻的目标状态矢量，构成初步估计的解集。
对于分组后的单组数据如图2所示，点E表示坐标系中心，定义观测点与目标之间的距离大小为ρi，则有：
Figure BDA0001999469440000122
其中，ti时刻卫星观测位置向量和目标对应的位置向量分别为
Figure BDA0001999469440000123
和
Figure BDA0001999469440000124
观测时刻t1、t2和t3对应的单位视线向量为
Figure BDA0001999469440000125
和
Figure BDA0001999469440000126
上式中有6个未知量：
Figure BDA0001999469440000127
和ρ1、ρ2、ρ3；r2为
Figure BDA0001999469440000128
的模长，所得目标位置如图3所示，目标速度如图4所示，其中，*号标记表示计算获得的位置解，+号标记表示目标实际位置解，▲标记表示观测卫星位置。根据Gauss方法，可求得
Figure BDA0001999469440000129
的初值以及
Figure BDA00019994694400001210
本实例中所得数据见表3，为实施例1中Gauss法求取目标位置、速度解群数据。
表3
Figure BDA00019994694400001211
Figure BDA0001999469440000131
S201.剔除初步估计的解集中不合理的解。
剔除原则为：剔除满足
Figure BDA0001999469440000132
的位置，速度数据。本实例中获取位置、速度解见表4，为实施例1中平滑前的目标位置误差数据。
表4
Figure BDA0001999469440000133
可表示如下：
Figure BDA0001999469440000141
Figure BDA0001999469440000142
其中，N表示经过剔除不合理数据后，剩余的位置或速度解的个数。
S202.将剩余的状态矢量解集拆分为位置分矢量解集和速度分矢量解集，分别进行分群，使得正确的解尽可能在同一个群中，得到位置分矢量解群和速度分矢量解群。
对所求的
Figure BDA0001999469440000143
分别进行k均值聚类，优选的，若所有时刻对应的位置解个数为1，则不需要进行k均值聚类，若其中任意时刻对应的位置解个数最大为2，则进行一次k为2的k均值聚类；若其中任意时刻对应的位置解个数最大为3，则进行1次k为3的k均值聚类，速度解聚类同上。
以位置解为例，聚类结果如图5所示，其中，*号标记点为类1，o号标记点为类2，+号标记点为目标实际位置，▲标记点为观测卫星位置，本实例中数据可见表5，为实施例1中平滑后的目标位置误差数据。
表5
Figure BDA0001999469440000144
Figure BDA0001999469440000151
实施例1采用Chauvenet's-criterion(肖维奈准则)判别进行异常数据消除
取置信概率为
Figure BDA0001999469440000152
时，Chauvenet's-criterion的置信概率公式可表示为：wn＝1+0.4In(n)。
对于待检测值xi，计算其与样本均值之间做差的绝对值，如若满足下式，则剔除当前待检测值：
Figure BDA0001999469440000153
其中，xi为待检测值，此处选取位置、速度向量在(x,y,z)三个方向上的分量分别进行三次检测，若任意一次满足则剔除，
Figure BDA0001999469440000156
为此类样本均值，wn为肖维勒准则的置信概率，Sx为此类样本标准差。
S203.对于位置分矢量解群和速度分矢量解群分别进行降噪预处理。
本发明优选基于Savitzky-Golay滤波方法进行降噪，处理如下：
Figure BDA0001999469440000154
使用上述五点三阶模板平滑后能够有效地在(x,y,z)三方向上分别降噪，平滑后的位置数据在解空间内将更加集中，以位置解群为例，平滑前后的结果如图6所示，其中，*标记点为平滑前目标位置，o标记点为平滑后目标位置。
S301.对降噪后的位置分矢量解群和速度分矢量解群进行拟合，构建各个时刻的状态矢量组合。
对于平滑降噪后的数据，在x,y,z三个方向上分别作二次多项式拟合，计算各个中间时刻
Figure BDA0001999469440000155
对应的目标位置在x,y,z分量上的估计值后可组合得到对应的位置修正估计值
Figure BDA0001999469440000161
与速度修正估计值
Figure BDA0001999469440000162
S302.根据各个时刻对应的状态矢量组合，生成目标轨道三维轨迹解集。
S303.根据观测平台测量状态，将目标轨道三维轨迹解集投影到瞬时观测像面得到二维轨迹解集。
寻找相同时间下对应的目标拟合位置进行组合，得到2×2×34对位置、速度时间组合
Figure BDA0001999469440000163
计算其对应的轨道根数
Figure BDA0001999469440000164
对应的目标轨道
Figure BDA0001999469440000165
以及目标轨道投影到像平面上的轨迹：
Figure BDA0001999469440000166
所得目标预测轨道如图7所示，其中，·号标记点表示估计轨道，o号标记点表示目标仿真轨道。
S401.对每条二维轨迹,计算导数误差与位置误差。
第m条二维轨迹的导数误差的计算公式如下：
Figure BDA0001999469440000167
Figure BDA0001999469440000168
Figure BDA0001999469440000169
Figure BDA00019994694400001610
Figure BDA00019994694400001611
其中，U″*，V″*，U'*，V'*为二维轨迹
Figure BDA00019994694400001612
的二次多项式拟合参数，
Figure BDA00019994694400001613
为估计轨迹l′m在tn时刻对应的像面坐标，m表示第m条二维轨迹，k为观测时刻总数。
相机参数(朝向)保持不变，可以得到各个时刻对目标观测的相面坐标
Figure BDA0001999469440000171
对于这些观测数据在U,V方向分别做二次多项式拟合，可以求得第二项、第四项。
第m条二维轨迹的位置误差的计算公式如下：
Figure BDA0001999469440000172
Figure BDA0001999469440000173
Figure BDA0001999469440000174
其中，
Figure BDA0001999469440000175
为预测值与观测值之间的像面距离，
Figure BDA0001999469440000176
为估计轨迹与观测轨迹交点和观测点之间的像面距离，tn为观测时刻，
Figure BDA0001999469440000177
为估计轨迹l′m在tn时刻对应的像面坐标，
Figure BDA0001999469440000178
为观测轨迹l*在tn时刻对应的像面坐标；
Figure BDA0001999469440000179
为估计轨迹l′m与观测轨迹交点的像面坐标。
S402.综合考虑导数误差与位置误差，选择最优二维轨迹，将其对应的轨道根数作为最佳估计值。
对于第m条轨道，分别对导数误差和位置误差从小到大进行排序，得到排名序号Rankvel,m和Rankdis,m，两者在1～Nsv之间，并计算两者之和作为最终误差评价：
Rankm＝Rankvel,m+Rankdis,m
对于最终误差评价，取其最小者mopt作为最佳预测轨迹：
Figure BDA00019994694400001710
其中，Nsv为各个时刻对应的目标位置、速度数据组合所求得的轨道根数的个数。
S403.输出最优二维轨迹对应轨道根数，完成初轨确定。
如图8所示，其中，·号标记点为所有估计轨道，▲表示优选的估计轨道，o号标记点为目标仿真轨道。
以上，仅为本申请较佳的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应该以权利要求的保护范围为准。

Claims (9)

1.一种基于Gauss解群优选的短弧初轨确定方法，其特征在于，该方法包括以下步骤：
S1.将观测数据进行分组，对每组数据，采用Gauss方法求出对应时刻的目标状态矢量，构成初步估计的解集；
S2.将初步估计的解集拆分为位置分矢量解集和速度分矢量解集，分别进行分群，得到位置分矢量解群和速度分矢量解群；
S3.基于位置分矢量解群和速度分矢量解群，生成估计的二维轨迹解集；
S4.采用轨迹优选的方法对每条估计的二维轨迹进行评估，计算最优估计的二维轨迹对应轨道根数，完成初轨确定；
步骤S4包括以下子步骤：
S401.对每条估计的二维轨迹，计算导数误差与位置误差；
第m条估计的二维轨迹的导数误差的计算公式如下：
Figure FDA0002590602010000011
Figure FDA0002590602010000012
Figure FDA0002590602010000013
Figure FDA0002590602010000014
Figure FDA0002590602010000015
第m条估计的二维轨迹的位置误差的计算公式如下：
Figure FDA0002590602010000021
Figure FDA0002590602010000022
Figure FDA0002590602010000023
S402.综合考虑导数误差与位置误差，选择最优估计的二维轨迹，将其对应的轨道根数作为最佳估计值；
S403.输出最优估计的二维轨迹对应轨道根数，完成初轨确定；
其中，
Figure FDA0002590602010000024
为估计的二维轨迹l′m在tn时刻对应的像面坐标，U″*,U′*,
Figure FDA00025906020100000211
V″*,V′*,
Figure FDA0002590602010000025
为观测的二维轨迹
Figure FDA0002590602010000026
的二次多项式拟合参数，
Figure FDA0002590602010000027
为观测的二维轨迹l*在tn时刻对应的像面坐标，m表示第m条估计的二维轨迹，k为观测时刻总数；
Figure FDA0002590602010000028
为预测值与观测值之间的像面距离，
Figure FDA0002590602010000029
为估计的二维轨迹与观测的二维轨迹交点和观测点之间的像面距离，tn为观测时刻；
Figure FDA00025906020100000210
为估计的二维轨迹l′m与观测的二维轨迹交点的像面坐标。
2.如权利要求1所述的短弧初轨确定方法，其特征在于，步骤S1包括以下子步骤：
S101.对观测数据进行分组，每3个时刻对应的矢量数据分为一组；
S102.对每组矢量数据，采用传统Gauss法求出对应时刻的目标状态矢量，构成初步估计的解集。
3.如权利要求2所述的短弧初轨确定方法，其特征在于，假定观测时刻为[t0,tk]，k+1≥3，对于[t0,tk]共k+1个时刻的输入数据，取间隔参数
Figure FDA0002590602010000031
其中，
Figure FDA0002590602010000032
表示向下取整，rate为间隔率，分组后的观测数据表示为
Figure FDA0002590602010000033
其中，
Figure FDA0002590602010000034
表示卫星观测位置矢量，
Figure FDA0002590602010000035
表示观测角度矢量，
Figure FDA0002590602010000036
Figure FDA0002590602010000037
4.如权利要求1所述的短弧初轨确定方法，其特征在于，步骤S2包括以下子步骤：
S201.剔除初步估计的解集中不合理的解；
S202.将剩余的状态矢量解集拆分为位置分矢量解集和速度分矢量解集，分别进行分群，使得正确的解尽可能在同一个群中，得到位置分矢量解群和速度分矢量解群；
S203.对于位置分矢量解群和速度分矢量解群分别进行降噪预处理。
5.如权利要求4所述的短弧初轨确定方法，其特征在于，剔除满足
Figure FDA0002590602010000038
或者
Figure FDA0002590602010000039
或者
Figure FDA00025906020100000310
任一条件的位置数据和速度数据，其中，
Figure FDA00025906020100000311
表示ts时刻的目标位置矢量，
Figure FDA00025906020100000312
表示ts时刻的目标速度矢量。
6.如权利要求4所述的短弧初轨确定方法，其特征在于，步骤S202中采用k-means聚类方法进行聚类分群，并采用Chauvenet's-criterion判别进行异常数据消除。
7.如权利要求1所述的短弧初轨确定方法，其特征在于，步骤S3包括以下子步骤：
S301.对降噪后的位置分矢量解群和速度分矢量解群进行拟合，构建各个时刻的状态矢量组合；
S302.根据各个时刻对应的状态矢量组合，生成目标轨道估计的三维轨迹解集；
S303.根据观测平台测量状态，将目标轨道估计的三维轨迹解集投影到瞬时观测像面，得到估计的二维轨迹解集。
8.如权利要求7所述的短弧初轨确定方法，其特征在于，步骤S402具体为：
对于第m条轨道，分别对导数误差和位置误差从小到大进行排序，得到排名序号Rankvel,m和Rankdis,m，两者在1～Nsv之间，并计算两者之和作为最终误差评价：
Rankm＝Rankvel,m+Rankdis,m
对于最终误差评价，取其最小者mopt作为最佳轨道根数估计值：
Figure FDA0002590602010000041
其中，Nsv为各个时刻对应的目标位置、速度数据组合所求得的轨道根数的个数。
9.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1至8任一项所述的短弧初轨确定方法。
CN201910207587.6A 2019-03-19 2019-03-19 一种基于Gauss解群优选的短弧初轨确定方法 Active CN110017832B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201910207587.6A CN110017832B (zh) 2019-03-19 2019-03-19 一种基于Gauss解群优选的短弧初轨确定方法
US16/762,505 US11662209B2 (en) 2019-03-19 2019-09-23 Short arc initial orbit determining method based on gauss solution cluster
PCT/CN2019/107186 WO2020186719A1 (zh) 2019-03-19 2019-09-23 一种基于Gauss解群优选的短弧初轨确定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910207587.6A CN110017832B (zh) 2019-03-19 2019-03-19 一种基于Gauss解群优选的短弧初轨确定方法

Publications (2)

Publication Number Publication Date
CN110017832A CN110017832A (zh) 2019-07-16
CN110017832B true CN110017832B (zh) 2020-10-16

Family

ID=67189639

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910207587.6A Active CN110017832B (zh) 2019-03-19 2019-03-19 一种基于Gauss解群优选的短弧初轨确定方法

Country Status (3)

Country Link
US (1) US11662209B2 (zh)
CN (1) CN110017832B (zh)
WO (1) WO2020186719A1 (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110017832B (zh) * 2019-03-19 2020-10-16 华中科技大学 一种基于Gauss解群优选的短弧初轨确定方法
CN111551183B (zh) * 2020-06-09 2021-12-21 中国人民解放军63921部队 基于天基光学观测数据的geo目标多点择优短弧定轨方法
CN113297745B (zh) * 2021-05-28 2022-09-02 中国人民解放军63921部队 一种基于短弧拟合位置的双弧段轨道改进方法
CN114387332B (zh) * 2022-01-17 2022-11-08 江苏省特种设备安全监督检验研究院 一种管道测厚方法及装置
CN115017386B (zh) * 2022-08-08 2022-11-08 中国人民解放军63921部队 基于聚类的观测数据与目标库根数关联方法和装置
CN115659196B (zh) * 2022-12-13 2023-06-23 中国人民解放军国防科技大学 基于非线性偏差演化的天基光学观测短弧关联与聚类方法
CN116202535B (zh) * 2022-12-28 2024-01-19 北京理工大学 一种初值智能优选的航天器仅测角超短弧初轨确定方法

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US3290933A (en) * 1965-10-18 1966-12-13 Robert L Lillestrand Navigation systems
US7471211B2 (en) * 1999-03-03 2008-12-30 Yamcon, Inc. Celestial object location device
US6237876B1 (en) * 2000-07-28 2001-05-29 Space Systems/Loral, Inc. Methods for using satellite state vector prediction to provide three-axis satellite attitude control
US7339731B2 (en) * 2005-04-20 2008-03-04 Meade Instruments Corporation Self-aligning telescope
US8718323B2 (en) * 2010-10-19 2014-05-06 Raytheon Company Batch detection association for enhanced target descrimination in dense detection environments
CN102968552B (zh) * 2012-10-26 2016-01-13 郑州威科姆科技股份有限公司 一种卫星轨道数据预估与修正方法
CN103364836A (zh) * 2013-07-23 2013-10-23 中国西安卫星测控中心 月球探测器落点预报方法
CN103927289B (zh) * 2014-04-23 2017-06-27 上海微小卫星工程中心 一种依据天基卫星测角资料确定低轨目标卫星初始轨道的方法
CN104457705B (zh) * 2014-12-12 2016-10-26 北京理工大学 基于天基自主光学观测的深空目标天体初定轨方法
CN104794268B (zh) * 2015-04-09 2017-12-26 中国科学院国家天文台 一种利用空间密度分布生成空间物体轨道的方法
CN105607058B (zh) * 2015-12-24 2018-05-01 中国科学院电子学研究所 利用geosar相位定标信息进行短弧段轨道精化的方法
CN107153209B (zh) * 2017-07-06 2019-07-30 武汉大学 一种短弧段低轨导航卫星实时精密定轨方法
CN108761507B (zh) * 2018-05-21 2020-07-03 中国人民解放军战略支援部队信息工程大学 基于短弧定轨和预报的导航卫星轨道快速恢复方法
CN110017832B (zh) * 2019-03-19 2020-10-16 华中科技大学 一种基于Gauss解群优选的短弧初轨确定方法

Also Published As

Publication number Publication date
CN110017832A (zh) 2019-07-16
US20210199439A1 (en) 2021-07-01
WO2020186719A1 (zh) 2020-09-24
US11662209B2 (en) 2023-05-30

Similar Documents

Publication Publication Date Title
CN110017832B (zh) 一种基于Gauss解群优选的短弧初轨确定方法
CN106803271B (zh) 一种视觉导航无人机的摄像机标定方法及装置
US20170131716A1 (en) Methods and apparatus to autonomously navigate a vehicle by selecting sensors from which to obtain measurements for navigation
CN110007317B (zh) 一种选星优化的高级接收机自主完好性监测方法
US10371787B2 (en) Emitter geolocation using sorted observations
EP2578995B1 (en) Modified Kalman filter for generation of attitude error corrections
CN110986968B (zh) 三维重建中实时全局优化和错误回环判断的方法及装置
CN110081881A (zh) 一种基于无人机多传感器信息融合技术的着舰引导方法
CA2920519C (en) Angles-only initial orbit determination (iod)
CN112444246B (zh) 高精度的数字孪生场景中的激光融合定位方法
CN110412868B (zh) 一种使用星间光学图像的非合作航天器轨道确定方法
CN112113582A (zh) 时间同步处理方法、电子设备及存储介质
CN110032201A (zh) 一种基于卡尔曼滤波的imu机载视觉姿态融合的方法
CN111998855B (zh) 光学望远镜共视观测确定空间目标初轨的几何方法及系统
CN111538339A (zh) 船舶航迹控制方法及装置
CN108225307A (zh) 一种惯性测量信息辅助的星图匹配方法
CN109870716B (zh) 定位方法和定位装置及计算机可读存储介质
CN111024124A (zh) 一种多传感器信息融合的组合导航故障诊断方法
CN112083457B (zh) 一种神经网络优化的imm卫星定位导航方法
Pastor et al. Initial orbit determination methods for track-to-track association
CN112269174A (zh) 基于交互式多模型信息融合的助推滑翔飞行器目标估计方法及系统
CN115143955A (zh) 基于天文测角数据确定地球同步轨道带航天器初轨的方法
CN110111356B (zh) 运动旋转物体的转动轴方向与转动角速度解算方法
CN113587926B (zh) 一种航天器空间自主交会对接相对导航方法
JP7246427B2 (ja) 画像を処理するための方法と装置、電子機器、コンピュータ可読記憶媒体及びコンピュータープログラム

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