WO2020186719A1 - 一种基于Gauss解群优选的短弧初轨确定方法 - Google Patents
一种基于Gauss解群优选的短弧初轨确定方法 Download PDFInfo
- Publication number
- WO2020186719A1 WO2020186719A1 PCT/CN2019/107186 CN2019107186W WO2020186719A1 WO 2020186719 A1 WO2020186719 A1 WO 2020186719A1 CN 2019107186 W CN2019107186 W CN 2019107186W WO 2020186719 A1 WO2020186719 A1 WO 2020186719A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- vector
- trajectory
- group
- observation
- sub
- 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.)
- Ceased
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
Definitions
- the present invention belongs to the field of astrodynamics, and more specifically, relates to a short arc initial orbit determination method based on Gauss solution group optimization.
- the traditional method for determining the initial orbit is to calculate the orbit of the satellite or space debris through ground observation data. These methods are usually based on the equation of motion of the two-body problem. Due to the short observation time of the space-based observing system and less available data, the traditional initial orbit determination method is not suitable for space-based observation data; in addition, if the curvature of the trajectory to be detected is too small, the initial orbit may not be determined.
- the current primary orbit determination methods mainly include: the traditional Laplace method, Gauss method and related improved algorithms.
- the shortcomings of the traditional Gauss method are: too many roots are obtained after solving high-order polynomials and need to be filtered; the improved Gauss method can avoid the problem of multiple roots, but it needs to choose a suitable initial value, otherwise it will cause the algorithm to diverge or eventually obtain a trivial solution. Wrong solution.
- the purpose of the present invention is to solve the technical problem that the prior art based on the traditional Gauss method for short arc initial orbit determination method has large errors and is difficult to effectively select from multiple sets of solutions.
- embodiments of the present invention provide a short arc initial orbit determination method based on Gauss solution group optimization, and the method includes the following steps:
- step S1 includes the following sub-steps:
- the traditional Gauss method is used to obtain the target state vector at the corresponding time to form a preliminary estimated solution set.
- the interval parameter among them Indicates rounding down
- rate is the interval rate
- the grouped observation data is expressed as among them, Represents the satellite observation position vector, Represents the observation angle vector,
- step S2 includes the following sub-steps:
- step S202 the k-means clustering method is used for clustering and grouping, and the Chauvenet's-criterion discrimination is used to eliminate abnormal data.
- step S3 includes the following sub-steps:
- step S4 includes the following sub-steps:
- U′′ * , V′′ * , U′ * , V′ * are two-dimensional trajectories
- t n is the observation time
- Is the observation trajectory l * corresponding image plane coordinates at t n To estimate the image plane coordinates of the intersection of the trajectory l′ m and the observation trajectory;
- step S402 is specifically:
- N sv is the number of track elements obtained by combining the target position and velocity data corresponding to each time.
- an embodiment of the present invention provides a computer-readable storage medium having a computer program stored on the computer-readable storage medium, and when the computer program is executed by a processor, the short arc initial path described in the first aspect is realized. Determine the method.
- the observation data is grouped and the traditional Gauss method is invoked for each group of data to increase the utilization of the data.
- the number of short arc target orbits finally determined by combining multiple groups of data will be closer to the true value.
- the present invention calculates the number of track roots corresponding to all roots after removing unreasonable data, noise reduction and smoothing, etc., and uses the method of trajectory priority to evaluate, and finally selects the highest evaluation index
- the number of orbital elements is used as the estimated solution, which can effectively solve the multi-root problem in the traditional Gauss method.
- FIG. 1 is a flowchart of a method for determining the initial orbit of a short arc based on Gauss solution group optimization according to an embodiment of the present invention
- FIG. 2 is a schematic diagram of space-based observation provided by an embodiment of the present invention.
- FIG. 3 is a solution group diagram of the target position obtained by the Gauss method according to an embodiment of the present invention.
- FIG. 5 is a position solution group diagram after smoothing by Savitzky-Golay provided by an embodiment of the present invention.
- FIG. 6 is a three-dimensional diagram of the orbit corresponding to each position solution obtained by calculating the number of orbits provided by an embodiment of the present invention
- FIG. 7 is a two-dimensional diagram of a three-dimensional orbit projected on an image surface with a solution corresponding to each position according to an embodiment of the present invention
- FIG. 8 is a three-dimensional diagram of the optimal orbit and the actual observation orbit provided by an embodiment of the present invention.
- Solution group the set of solutions.
- the present invention provides a short arc initial orbit determination method based on Gauss solution group optimization, which includes the following steps:
- Step S1 Group the observation data, and for each group of data, the Gauss method is used to obtain the target state vector at the corresponding time to form a preliminary estimated solution set.
- Observation data is observation data of space-based monocular imaging platform or ground observation data.
- the observation data includes: the observation position corresponding to each observation time (the position of the platform itself relative to the geocentric coordinate system), the observation angle vector, and the variable Represents the satellite observation position vector at time t s , Indicates the observation angle vector at time t s .
- the position of the observation satellite, the angle of the observation camera, and the observation imaging data at the same time constitute a pair of data, and multiple time data constitute a set of data.
- the present invention first needs to group the observation data, and the vector data corresponding to every three moments are grouped into one group. Assuming that the observation time is [t 0 ,t k ](k+1 ⁇ 3), for the input data of [t 0 ,t k ] total k+1 time, take the interval parameter among them, Represents rounding down, rate is the interval rate, usually 0.2.
- the grouped observation data is expressed as among them, Represents the satellite observation position vector, Represents the observation angle vector,
- the traditional Gauss method is used to obtain the target state vector at the corresponding time to form a preliminary estimated solution set.
- the present invention considers using the Gauss method multiple times to increase the utilization of the data.
- the number of short-arc target tracks finally determined by combining multiple sets of data will be closer to reality value.
- the Gauss method solves the target position vector and target velocity vector at the intermediate time s2 by inputting the observation position vector and the observation angle vector corresponding to three different moments.
- Targets refer to unpowered targets, which can be debris or out-of-control satellites.
- every three of a series of moments are divided into a group, and each group of data is solved by a traditional Gauss method to obtain the target position vector and velocity vector solution at the middle moment of the group of data.
- the number of solutions obtained for each set of vector data may be 1, 2, or 3.
- the Gauss method is called multiple times to obtain the target vector corresponding to different moments to form a set of preliminary estimated solutions.
- Step S2 Split the preliminarily estimated solution set into a position sub-vector solution set and a velocity sub-vector solution set, and perform grouping respectively to obtain the position sub-vector solution group and the velocity sub-vector solution group.
- Eliminate satisfaction (Earth radius) or (Geosynchronous geostationary orbit radius) or (Second Universe Velocity) Position data and velocity data of any condition, in which, Represents the target position vector at time t s , Indicates the target velocity vector at time t s .
- a solution that satisfies any of the above three conditions is a solution that does not conform to common sense. Therefore, the target position data whose target position is less than the radius of the earth or the radius of the geosynchronous geostationary orbit and the target speed data whose target position is greater than the second cosmic velocity are eliminated.
- the present invention preferably uses the k-means clustering method to perform clustering (k is 2 or 3), and uses Chauvenet's-criterion to eliminate abnormal data.
- the abnormal data refers to the performance in the same group Unusual data, for example, outlier data.
- the present invention preferably performs noise reduction based on the Savitzky-Golay filtering method.
- Step S3. Generate a two-dimensional trajectory solution set based on the position sub-vector solution group and the velocity sub-vector solution group.
- the target position component vector fitting and target velocity component vector fitting are performed.
- the combination of velocity vector and position vector corresponding to time t s is
- the 3D trajectory solution set can be determined by the number of trajectories Among them, a is the semi-major axis, e is the first eccentricity, i is the orbital inclination, ⁇ is the right ascension of the ascending node, ⁇ is the argument of perigee, and M is the mean anomaly.
- a is the semi-major axis
- e is the first eccentricity
- i the orbital inclination
- ⁇ the right ascension of the ascending node
- ⁇ is the argument of perigee
- M is the mean anomaly.
- Step S4 Use the method of trajectory optimization to evaluate each two-dimensional trajectory, calculate the number of trajectories corresponding to the optimal two-dimensional trajectory, and complete the initial trajectory determination.
- U′′ * , V′′ * , U′ * , V′ * are two-dimensional trajectories
- the camera parameters (orientation) remain unchanged, and the phase plane coordinates of the target observation at each time can be obtained
- the second and fourth terms can be obtained by fitting quadratic polynomials in the U and V directions respectively.
- t n is the observation time
- Is the observed trajectory l * corresponding image plane coordinates at t n To estimate the image plane coordinates of the intersection point of the trajectory l' m and the observation trajectory.
- N sv is the number of track elements obtained by combining the target position and velocity data corresponding to each time.
- the invention can effectively solve multiple problems that are difficult to solve by the traditional Gauss method, make full use of the original observation data, and use the concept of group solution and the optimization method to improve the accuracy of initial orbit determination. Under the same input error condition, this technical solution It performs better than other solutions for short-term observation data, that is, the present invention is more suitable for short-arc observations than other inventions.
- Step S102 For each group of vector data, the traditional Gauss method is used to obtain the target state vector at the corresponding time to form a preliminary estimated solution set.
- point E represents the center of the coordinate system, defining the distance between the observation point and the target as ⁇ i , then:
- the satellite observation position vector and the position vector corresponding to the target at t i are respectively with
- the unit line of sight vector corresponding to the observation time t 1 , t 2 and t 3 is with
- r 2 is The target position is shown in Fig. 3, and the target speed is shown in Fig. 4, where the * mark indicates the calculated position solution, the + mark indicates the actual target position solution, and the ⁇ mark indicates the observation satellite position.
- Table 3 is the Gauss method in Example 1 to obtain the target position and velocity solution group data.
- the elimination principle is: eliminate satisfaction
- the location, speed data The position and velocity solutions obtained in this example are shown in Table 4, which is the target position error data before smoothing in Embodiment 1.
- N represents the number of remaining position or velocity solutions after removing unreasonable data.
- K-means clustering is performed separately.
- the clustering result is shown in Figure 5.
- the * mark is class 1
- the o mark is class 2
- the + mark is the actual position of the target
- the ⁇ mark is the observation satellite position.
- Table 5 is the target position error data after smoothing in Example 1.
- x i For the value x i to be detected, calculate the absolute value of the difference between it and the sample mean. If the following formula is satisfied, the current value to be detected is eliminated: Among them, x i is the value to be detected.
- the components of the position and velocity vector in the three directions (x, y, z) are selected for three detections respectively, and if any one is satisfied, it will be eliminated. Is the mean value of such a sample, w n is the confidence probability of the Shawville criterion, and S x is the standard deviation of such a sample.
- the present invention preferably performs noise reduction based on the Savitzky-Golay filtering method, and the processing is as follows:
- the obtained target predicted trajectory is shown in Fig. 7, where the marked point marked with ⁇ represents the estimated trajectory, and the marked point marked o represents the simulated target trajectory.
- U′′ * , V′′ * , U′ * , V′ * are two-dimensional trajectories
- the camera parameters (orientation) remain unchanged, and the phase plane coordinates of the target observation at each time can be obtained
- the second and fourth terms can be obtained by fitting quadratic polynomials in the U and V directions respectively.
- t n is the observation time
- Is the observation trajectory l * corresponding image plane coordinates at t n To estimate the image plane coordinates of the intersection point of the trajectory l' m and the observation trajectory.
- N sv is the number of track elements obtained by combining the target position and velocity data corresponding to each time.
- the mark points marked with ⁇ are all estimated orbits, ⁇ is the preferred estimated track, and marked with o is the target simulation track.
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Astronomy & Astrophysics (AREA)
- Aviation & Aerospace Engineering (AREA)
- Chemical & Material Sciences (AREA)
- Combustion & Propulsion (AREA)
- Image Analysis (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
一种基于Gauss解群优选的短弧初轨确定方法,属于天文动力学领域。包括:将观测数据进行分组,对每组数据,采用Gauss方法求出对应时刻的目标状态矢量,构成初步估计的解集(S1);将初步估计的解集拆分为位置分矢量解集和速度分矢量解集,分别进行分群,得到位置分矢量解群和速度分矢量解群(S2);基于位置分矢量解群和速度分矢量解群,生成二维轨迹解集(S3);采用轨迹优选的方法对每条二维轨迹进行评估,计算最优二维轨迹对应轨道根数,完成初轨确定(S4)。对观测数据进行分组,对每组数据调用Gauss,加大数据的利用程度,结合多组数据最终确定的短弧目标轨道根数更贴近真实。采用Gauss解群优选进行轨迹评估,最优二维轨迹的轨道根数作为估计解,解决了多根问题。
Description
本发明属于天文动力学领域,更具体地,涉及一种基于Gauss解群优选的短弧初轨确定方法。
随着天文动力学的发展,初轨确定技术已成为当今最重要的技术之一。传统的初轨确定方法是通过地面观测数据计算卫星或空间碎片的轨道,这些方法通常是基于两体问题的运动方程。由于天基观测系统观测时刻短、可获取的数据较少,导致传统的初轨确定方法不适用于天基观测数据;此外,如果待检测轨迹的曲率过小,则有可能无法确定初轨。
当前的初轨确定方法主要包括:传统的Laplace法、Gauss法及相关改进算法。传统Gauss法的缺陷是:求解高阶多项式后所得的根太多,需要进行筛选;改进Gauss法可以避免重根问题,但需要选取合适的初值,否则会使得算法迭代发散或最终取的平凡解甚至错误解。
当前研究资料表明,在基于地面观测数据的初轨确定领域,已有许多性能优异的方法,而在基于天基观测数据的初轨确定领域,相关研究算法的时间间隔要求都比较长,对于短弧初轨确定适用性不高,本方法可解决传统Gauss方法中存在的多根问题,并适用于短弧初轨确定任务。
【发明内容】
针对现有技术的缺陷,本发明的目的在于解决现有技术基于传统Gauss方法的短弧初轨确定方法误差较大、难以从多组解中有效选择的技术问题。
为实现上述目的,第一方面,本发明实施例提供了一种基于Gauss解 群优选的短弧初轨确定方法,该方法包括以下步骤:
S1.将观测数据进行分组,对每组数据,采用Gauss方法求出对应时刻的目标状态矢量,构成初步估计的解集;
S2.将初步估计的解集拆分为位置分矢量解集和速度分矢量解集,分别进行分群,得到位置分矢量解群和速度分矢量解群;
S3.基于位置分矢量解群和速度分矢量解群,生成二维轨迹解集;
S4.采用轨迹优选的方法对每条二维轨迹进行评估,计算最优二维轨迹对应轨道根数,完成初轨确定。
具体地,步骤S1包括以下子步骤:
S101.对观测数据进行分组,每3个时刻对应的矢量数据分为一组;
S102.对每组矢量数据,采用传统Gauss法求出对应时刻的目标状态矢量,构成初步估计的解集。
具体地,假定观测时刻为[t
0,t
k],k+1≥3,对于[t
0,t
k]共k+1个时刻的输入数据,取间隔参数
其中,
表示向下取整,rate为间隔率,分组后的观测数据表示为
其中,
表示卫星观测位置矢量,
表示观测角度矢量,
具体地,步骤S2包括以下子步骤:
S201.剔除初步估计的解集中不合理的解;
S202.将剩余的状态矢量解集拆分为位置分矢量解集和速度分矢量解集,分别进行分群,使得正确的解尽可能在同一个群中,得到位置分矢量解群和速度分矢量解群;
S203.对于位置分矢量解群和速度分矢量解群分别进行降噪预处理。
具体地,步骤S202中采用k-means聚类方法进行聚类分群,并采用Chauvenet′s-criterion判别进行异常数据消除。
具体地,步骤S3包括以下子步骤:
S301.对降噪后的位置分矢量解群和速度分矢量解群进行拟合,构建各个时刻的状态矢量组合;
S302.根据各个时刻对应的状态矢量组合,生成目标轨道三维轨迹解集;
S303.根据观测平台测量状态,将目标轨道三维轨迹解集投影到瞬时观测像面,得到二维轨迹解集。
具体地,步骤S4包括以下子步骤:
S401.对每条二维轨迹,计算导数误差与位置误差;
第m条二维轨迹的导数误差的计算公式如下:
第m条二维轨迹的位置误差的计算公式如下:
其中,
为预测值与观测值之间的像面距离,
为估计轨迹与观测轨迹交点和观测点之间的像面距离,t
n为观测时刻,
为估计轨迹l′
m在t
n时刻对应的像面坐标,
为观测轨迹l
*在t
n时刻对应的像面坐标;
为估计轨迹l′
m与观测轨迹交点的像面坐标;
S402.综合考虑导数误差与位置误差,选择最优二维轨迹,将其对应的轨道根数作为最佳估计值;
S403.输出最优二维轨迹对应轨道根数,完成初轨确定。
具体地,步骤S402具体为:
对于第m条轨道,分别对导数误差和位置误差从小到大进行排序,得到排名序号Rank
vel,m和Rank
dis,m,两者在1~N
sv之间,并计算两者之和作为最终误差评价:
Rank
m=Rank
vel,m+Rank
dis,m
对于最终误差评价,取其最小者m
opt作为最佳预测轨迹:
其中,N
sv为各个时刻对应的目标位置、速度数据组合所求得的轨道根数的个数。
第二方面,本发明实施例提供了一种计算机可读存储介质,该计算机可读存储介质上存储有计算机程序,该计算机程序被处理器执行时实现上述第一方面所述的短弧初轨确定方法。
总体而言,通过本发明所构思的以上技术方案与现有技术相比,具有以下有益效果:
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个时刻对应的矢量数据分为一组。
观测数据是天基单目成像平台观测数据或者地面观测数据。所述观测数据包括:各个观测时刻对应的观测位置(平台本身相对地心坐标系位置)、观测角度矢量,分别用变量
表示t
s时刻卫星观测位置矢量、
表示t
s时刻观测角度矢量。同一时刻下观测卫星的位置、观测相机的角度以及观测成像数据共同构成一对数据,多个时刻数据构成一组数据。
为了后续构建解群,本发明首先需要对观测数据进行分组,每3个时刻对应的矢量数据分为一组。假定观测时刻为[t
0,t
k](k+1≥3),对于[t
0,t
k]共k+1个时刻的输入数据,取间隔参数
其中,
表示向下取整,rate为间隔率,通常取值0.2。分组后的观测数据表示为
其中,
表示卫星观测位置矢量,
表示观测角度矢量,
S102.对每组矢量数据,采用传统Gauss法求出对应时刻的目标状态矢量,构成初步估计的解集。
由于传统Gauss方法输入仅为3个数据,所输入的信息有限,本发明考虑多次使用Gauss方法以加大数据的利用程度,结合多组数据最终确定的短弧目标轨道根数会更加贴近真实值。
Gauss法是通过输入三个不同时刻对应的观测位置向量和观测角度向量,求解中间时刻s2的目标位置矢量与目标速度矢量。目标指无动力目标,可以是碎片或者失控卫星等。按照本发明方法将一系列时刻中每三个分为一组,每组数据均采用一次传统Gauss法求解得到该组数据中间时刻的目标位置矢量和速度矢量解。每组矢量数据得到的解的个数可能为1、2、3。多次调用Gauss方法,求出不同时刻对应的目标矢量构成初步估计的解的集合。
步骤S2.将初步估计的解集拆分为位置分矢量解集和速度分矢量解集,分别进行分群,得到位置分矢量解群和速度分矢量解群。
S201.剔除初步估计的解集中不合理的解。
很显然,满足上述三个条件任一个的解,是不符合常理的解。因此,剔除所求目标位置距离地心小于地球半径或大于地球同步静止轨道半径的目标位置数据以及大于第二宇宙速度的目标速度数据。
S202.将剩余的状态矢量解集拆分为位置分矢量解集和速度分矢量解集,分别进行分群,使得正确的解尽可能在同一个群中,得到位置分矢量解群和速度分矢量解群。
为了将正确的根和错误的根分开,尽量使得正确的根在同一个群,对位置分矢量解集和速度分矢量解集分别进行分群。本发明优选k-means聚类方法进行聚类分群(k取值2或者3),并采用Chauvenet′s-criterion(肖维奈准则)判别进行异常数据消除,异常数据是指同一个群中表现异常的数据,例如,离群数据。
若每个时刻对应的位置解个数小于等于1,则不需要进行k均值聚类;若其中任意时刻对应的位置解个数最大为2,则进行一次k=2的k均值聚类;若其中任意时刻对应的位置解个数最大为3,则进行一次k=3的k均值聚类。速度分矢量解集聚类同上。
S203.对于位置分矢量解群和速度分矢量解群分别进行降噪预处理。
本发明优选基于Savitzky-Golay滤波方法进行降噪。
所有根经过剔除不合理数据、降噪平滑等处理方式后,能够提升轨迹优选的精确度。
步骤S3.基于位置分矢量解群和速度分矢量解群,生成二维轨迹解集。
S301.对降噪后的位置分矢量解群和速度分矢量解群进行拟合,构建各个时刻的状态矢量组合。
S302.根据各个时刻对应的状态矢量组合,生成目标轨道三维轨迹解集。
通过某一时刻t
s对应的速度矢量和位置矢量组合
计算出
对应的轨道根数
通过轨道根数可以确定三维轨迹解集
其中,a为半长轴,e为第一偏心率,i为轨道倾角,Ω为升交点赤经,ω为近地点幅角,M为平近点角,
为观测轨迹
在t
n时刻对应的三维坐标,m=1,2,...,N
sv, n=0,1,...,k,N
sv为各个时刻对应的目标位置、速度数据组合所求得的轨道根数的总个数。
S303.根据观测平台测量状态,将目标轨道三维轨迹解集投影到瞬时观测像面,得到二维轨迹解集。
目标轨道三维轨迹解集
投影到瞬时观测像面上,即可得到二维轨迹解集
其中,
为第m条三维观测轨迹,
为观测轨迹
在在t
n时刻对应的三维坐标,l′
m为第m条估计轨迹,
为估计轨迹l′
m在t
n时刻对应的像面坐标,n=0,1,...,k,m=1,2,...,N
sv。
步骤S4.采用轨迹优选的方法对每条二维轨迹进行评估,计算最优二维轨迹对应轨道根数,完成初轨确定。
S401.对每条二维轨迹,计算导数误差与位置误差。
第m条二维轨迹的导数误差的计算公式如下:
第m条二维轨迹的位置误差的计算公式如下:
其中,
为预测值与观测值之间的像面距离,
为估计轨迹与观测轨迹交点和观测点之间的像面距离,t
n为观测时刻,
为估计轨迹l′
m在t
n时刻对应的像面坐标,
为观测轨迹l
*在t
n时刻对应的像面坐标,
为估计轨迹l′
m与观测轨迹交点的像面坐标。
S402.综合考虑导数误差与位置误差,选择最优二维轨迹,将其对应的轨道根数作为最佳估计值。
对于第m条轨道,分别对导数误差和位置误差从小到大进行排序,得到排名序号Rank
vel,m和Rank
dis,m,两者在1~N
sv之间,并计算两者之和作为最终误差评价:
Rank
m=Rank
vel,m+Rank
dis,m
对于最终误差评价,取其最小者m
opt作为最优二维轨迹:
其中,N
sv为各个时刻对应的目标位置、速度数据组合所求得的轨道根数的个数。
S403.输出最优二维轨迹对应轨道根数,完成初轨确定。
适用于单一天基成像观测平台,短时间观测条件(大于5s)下的远距离无动力目标,如碎片、失控卫星等的初轨确定。本发明能够有效解决传统Gauss方法难以解决的多根问题,对于原始观测数据的利用更加充分, 利用解群的概念和优选方法提高了初轨定轨精度,在相同输入误差条件下,本技术方案对短时间观测数据比其他方案表现更佳,即本发明较其他发明对于短弧观测更适用。
实施例1
如表1所示,为实施例1中观测卫星与目标轨道根数数据。
表1
S101.对观测数据进行分组,每3个时刻对应的矢量数据分为一组。
对于[t
0,t
170]共171个时刻的输入数据,取间隔参数Nrate=34,将观测数据
分为一组,其中,
表示向下取整,rate为间隔率取rate=0.2,
表示卫星观测位置向量,
表示观测角度向量。
对于表1数据下的观测卫星与目标轨道根数,仿真得到表2中观测卫星与目标部分位置数据。
表2
步骤S102.对每组矢量数据,采用传统Gauss法求出对应时刻的目标状态矢量,构成初步估计的解集。
对于分组后的单组数据如图2所示,点E表示坐标系中心,定义观测点与目标之间的距离大小为ρ
i,则有:
其中,t
i时刻卫星观测位置向量和目标对应的位置向量分别为
和
观测时刻t
1、t
2和t
3对应的单位视线向量为
和
上式中有6个未知量:
和ρ
1、ρ
2、ρ
3;r
2为
的模长,所得目标位置如图3所示,目标速度如图4所示,其中,*号标记表示计算获得的位置解,+号标记表示目标实际位置解,▲标记表示观测卫星位置。根据Gauss方法,可求得
的初值以及
本实例中所得数据见表3,为实施例1中Gauss法求取目标位置、速度解群数据。
表3
S201.剔除初步估计的解集中不合理的解。
表4
可表示如下:
其中,N表示经过剔除不合理数据后,剩余的位置或速度解的个数。
S202.将剩余的状态矢量解集拆分为位置分矢量解集和速度分矢量解集,分别进行分群,使得正确的解尽可能在同一个群中,得到位置分矢量解群和速度分矢量解群。
对所求的
分别进行k均值聚类,优选的,若所有时刻对应的位置解个数为1,则不需要进行k均值聚类,若其中任意时刻对应的位置解个数最大为2,则进行一次k为2的k均值聚类;若其中任意时刻对应的位置解个数最大为3,则进行1次k为3的k均值聚类,速度解聚类同上。
以位置解为例,聚类结果如图5所示,其中,*号标记点为类1,o号标记点为类2,+号标记点为目标实际位置,▲标记点为观测卫星位置,本实例中数据可见表5,为实施例1中平滑后的目标位置误差数据。
表5
实施例1采用Chauvenet′s-criterion(肖维奈准则)判别进行异常数据消除
对于待检测值x
i,计算其与样本均值之间做差的绝对值,如若满足下式,则剔除当前待检测值:
其中,x
i为待检测值,此处选取位置、速度向量在(x,y,z)三个方向上的分量分别进行三次检测,若任意一次满足则剔除,
为此类样本均值,w
n为肖维勒准则的置信概率,S
x为此类样本标准差。
S203.对于位置分矢量解群和速度分矢量解群分别进行降噪预处理。
本发明优选基于Savitzky-Golay滤波方法进行降噪,处理如下:
使用上述五点三阶模板平滑后能够有效地在(x,y,z)三方向上分别降噪,平滑后的位置数据在解空间内将更加集中,以位置解群为例,平滑前后的结果如图6所示,其中,*标记点为平滑前目标位置,o标记点为平滑后目标位置。
S301.对降噪后的位置分矢量解群和速度分矢量解群进行拟合,构建各个时刻的状态矢量组合。
S302.根据各个时刻对应的状态矢量组合,生成目标轨道三维轨迹解集。
S303.根据观测平台测量状态,将目标轨道三维轨迹解集投影到瞬时观测像面得到二维轨迹解集。
寻找相同时间下对应的目标拟合位置进行组合,得到2×2×34对位置、速度时间组合
计算其对应的轨道根数
对应的目标轨道
以及目标轨道投影到像平面上的轨迹:
所得目标预测轨道如图7所示,其中,·号标记点表示估计轨道,o号标记点表示目标仿真轨道。
S401.对每条二维轨迹,计算导数误差与位置误差。
第m条二维轨迹的导数误差的计算公式如下:
第m条二维轨迹的位置误差的计算公式如下:
其中,
为预测值与观测值之间的像面距离,
为估计轨迹与观测轨迹交点和观测点之间的像面距离,t
n为观测时刻,
为估计轨迹l′
m在t
n时刻对应的像面坐标,
为观测轨迹l
*在t
n时刻对应的像面坐标;
为估计轨迹l′
m与观测轨迹交点的像面坐标。
S402.综合考虑导数误差与位置误差,选择最优二维轨迹,将其对应的轨道根数作为最佳估计值。
对于第m条轨道,分别对导数误差和位置误差从小到大进行排序,得到排名序号Rank
vel,m和Rank
dis,m,两者在1~N
sv之间,并计算两者之和作为最终误差评价:
Rank
m=Rank
vel,m+Rank
dis,m
对于最终误差评价,取其最小者m
opt作为最佳预测轨迹:
其中,N
sv为各个时刻对应的目标位置、速度数据组合所求得的轨道根数的个数。
S403.输出最优二维轨迹对应轨道根数,完成初轨确定。
如图8所示,其中,·号标记点为所有估计轨道,▲表示优选的估计轨道,o号标记点为目标仿真轨道。
以上,仅为本申请较佳的具体实施方式,但本申请的保护范围并不局 限于此,任何熟悉本技术领域的技术人员在本申请揭露的技术范围内,可轻易想到的变化或替换,都应涵盖在本申请的保护范围之内。因此,本申请的保护范围应该以权利要求的保护范围为准。
Claims (10)
- 一种基于Gauss解群优选的短弧初轨确定方法,其特征在于,该方法包括以下步骤:S1.将观测数据进行分组,对每组数据,采用Gauss方法求出对应时刻的目标状态矢量,构成初步估计的解集;S2.将初步估计的解集拆分为位置分矢量解集和速度分矢量解集,分别进行分群,得到位置分矢量解群和速度分矢量解群;S3.基于位置分矢量解群和速度分矢量解群,生成二维轨迹解集;S4.采用轨迹优选的方法对每条二维轨迹进行评估,计算最优二维轨迹对应轨道根数,完成初轨确定。
- 如权利要求1所述的短弧初轨确定方法,其特征在于,步骤S1包括以下子步骤:S101.对观测数据进行分组,每3个时刻对应的矢量数据分为一组;S102.对每组矢量数据,采用传统Gauss法求出对应时刻的目标状态矢量,构成初步估计的解集。
- 如权利要求1所述的短弧初轨确定方法,其特征在于,步骤S2包括以下子步骤:S201.剔除初步估计的解集中不合理的解;S202.将剩余的状态矢量解集拆分为位置分矢量解集和速度分矢量解集,分别进行分群,使得正确的解尽可能在同一个群中,得到位置分矢量解群和速度分矢量解群;S203.对于位置分矢量解群和速度分矢量解群分别进行降噪预处理。
- 如权利要求4所述的短弧初轨确定方法,其特征在于,步骤S202中采用k-means聚类方法进行聚类分群,并采用Chauvenet′s-criterion判别进行异常数据消除。
- 如权利要求1所述的短弧初轨确定方法,其特征在于,步骤S3包括以下子步骤:S301.对降噪后的位置分矢量解群和速度分矢量解群进行拟合,构建各个时刻的状态矢量组合;S302.根据各个时刻对应的状态矢量组合,生成目标轨道三维轨迹解集;S303.根据观测平台测量状态,将目标轨道三维轨迹解集投影到瞬时观测像面,得到二维轨迹解集。
- 如权利要求1所述的短弧初轨确定方法,其特征在于,步骤S4包括以下子步骤:S401.对每条二维轨迹,计算导数误差与位置误差;第m条二维轨迹的导数误差的计算公式如下:第m条二维轨迹的位置误差的计算公式如下:其中, 为预测值与观测值之间的像面距离, 为估计轨迹与观测轨迹交点和观测点之间的像面距离,t n为观测时刻, 为估计轨迹l′ m在t n时刻对应的像面坐标, 为观测轨迹l *在t n时刻对应的像面坐标; 为估计轨迹l′ m与观测轨迹交点的像面坐标;S402.综合考虑导数误差与位置误差,选择最优二维轨迹,将其对应的轨道根数作为最佳估计值;S403.输出最优二维轨迹对应轨道根数,完成初轨确定。
- 一种计算机可读存储介质,其特征在于,所述计算机可读存储介质上存储有计算机程序,所述计算机程序被处理器执行时实现如权利要求1至9任一项所述的短弧初轨确定方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| US16/762,505 US11662209B2 (en) | 2019-03-19 | 2019-09-23 | Short arc initial orbit determining method based on gauss solution cluster |
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910207587.6A CN110017832B (zh) | 2019-03-19 | 2019-03-19 | 一种基于Gauss解群优选的短弧初轨确定方法 |
| CN201910207587.6 | 2019-03-19 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| WO2020186719A1 true WO2020186719A1 (zh) | 2020-09-24 |
Family
ID=67189639
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| PCT/CN2019/107186 Ceased WO2020186719A1 (zh) | 2019-03-19 | 2019-09-23 | 一种基于Gauss解群优选的短弧初轨确定方法 |
Country Status (3)
| Country | Link |
|---|---|
| US (1) | US11662209B2 (zh) |
| CN (1) | CN110017832B (zh) |
| WO (1) | WO2020186719A1 (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113297745A (zh) * | 2021-05-28 | 2021-08-24 | 中国人民解放军63921部队 | 一种基于短弧拟合位置的双弧段轨道改进方法 |
| CN115828037A (zh) * | 2022-11-24 | 2023-03-21 | 上海卫星工程研究所 | 天基光学监测平台空间碎片初始轨道确定方法及系统 |
Families Citing this family (12)
| 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目标多点择优短弧定轨方法 |
| CN114387332B (zh) * | 2022-01-17 | 2022-11-08 | 江苏省特种设备安全监督检验研究院 | 一种管道测厚方法及装置 |
| CN115017386B (zh) * | 2022-08-08 | 2022-11-08 | 中国人民解放军63921部队 | 基于聚类的观测数据与目标库根数关联方法和装置 |
| CN115837992B (zh) * | 2022-11-25 | 2025-07-11 | 上海卫星工程研究所 | 面向空间碎片的天基光学观测初轨关联方法和系统 |
| CN115659196B (zh) * | 2022-12-13 | 2023-06-23 | 中国人民解放军国防科技大学 | 基于非线性偏差演化的天基光学观测短弧关联与聚类方法 |
| CN116202535B (zh) * | 2022-12-28 | 2024-01-19 | 北京理工大学 | 一种初值智能优选的航天器仅测角超短弧初轨确定方法 |
| CN117310765B (zh) * | 2023-07-07 | 2026-02-27 | 哈尔滨工业大学 | 基于多星观测的单目标轨道确定方法、设备、介质和产品 |
| CN118520674B (zh) * | 2024-05-15 | 2024-10-01 | 中国人民解放军军事航天部队航天工程大学 | 基于最小容许域的光学短弧初定轨的概率密度分布分析方法、系统、设备和介质 |
| CN118797450B (zh) * | 2024-09-14 | 2024-12-20 | 中科星图测控技术股份有限公司 | 一种基于天基光学测量数据处理与识别方法 |
| CN118936482B (zh) * | 2024-09-26 | 2025-10-17 | 哈尔滨工业大学 | 一种基于粒子群算法的拉普拉斯极短弧定轨算法 |
| CN120573285B (zh) * | 2025-07-31 | 2025-10-21 | 中国星网网络创新研究院有限公司 | 一种低轨目标定轨方法及装置 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
| CN103364836A (zh) * | 2013-07-23 | 2013-10-23 | 中国西安卫星测控中心 | 月球探测器落点预报方法 |
| CN103927289A (zh) * | 2014-04-23 | 2014-07-16 | 上海微小卫星工程中心 | 一种依据天基卫星测角资料确定低轨目标卫星初始轨道的方法 |
| CN104457705A (zh) * | 2014-12-12 | 2015-03-25 | 北京理工大学 | 基于天基自主光学观测的深空目标天体初定轨方法 |
| CN107153209A (zh) * | 2017-07-06 | 2017-09-12 | 武汉大学 | 一种短弧段低轨导航卫星实时精密定轨方法 |
| CN110017832A (zh) * | 2019-03-19 | 2019-07-16 | 华中科技大学 | 一种基于Gauss解群优选的短弧初轨确定方法 |
Family Cites Families (8)
| 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 |
| 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 | 郑州威科姆科技股份有限公司 | 一种卫星轨道数据预估与修正方法 |
| CN104794268B (zh) * | 2015-04-09 | 2017-12-26 | 中国科学院国家天文台 | 一种利用空间密度分布生成空间物体轨道的方法 |
| CN105607058B (zh) * | 2015-12-24 | 2018-05-01 | 中国科学院电子学研究所 | 利用geosar相位定标信息进行短弧段轨道精化的方法 |
| CN108761507B (zh) * | 2018-05-21 | 2020-07-03 | 中国人民解放军战略支援部队信息工程大学 | 基于短弧定轨和预报的导航卫星轨道快速恢复方法 |
-
2019
- 2019-03-19 CN CN201910207587.6A patent/CN110017832B/zh active Active
- 2019-09-23 US US16/762,505 patent/US11662209B2/en active Active
- 2019-09-23 WO PCT/CN2019/107186 patent/WO2020186719A1/zh not_active Ceased
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
| CN103364836A (zh) * | 2013-07-23 | 2013-10-23 | 中国西安卫星测控中心 | 月球探测器落点预报方法 |
| CN103927289A (zh) * | 2014-04-23 | 2014-07-16 | 上海微小卫星工程中心 | 一种依据天基卫星测角资料确定低轨目标卫星初始轨道的方法 |
| CN104457705A (zh) * | 2014-12-12 | 2015-03-25 | 北京理工大学 | 基于天基自主光学观测的深空目标天体初定轨方法 |
| CN107153209A (zh) * | 2017-07-06 | 2017-09-12 | 武汉大学 | 一种短弧段低轨导航卫星实时精密定轨方法 |
| CN110017832A (zh) * | 2019-03-19 | 2019-07-16 | 华中科技大学 | 一种基于Gauss解群优选的短弧初轨确定方法 |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113297745A (zh) * | 2021-05-28 | 2021-08-24 | 中国人民解放军63921部队 | 一种基于短弧拟合位置的双弧段轨道改进方法 |
| CN115828037A (zh) * | 2022-11-24 | 2023-03-21 | 上海卫星工程研究所 | 天基光学监测平台空间碎片初始轨道确定方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| US20210199439A1 (en) | 2021-07-01 |
| CN110017832B (zh) | 2020-10-16 |
| US11662209B2 (en) | 2023-05-30 |
| CN110017832A (zh) | 2019-07-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| WO2020186719A1 (zh) | 一种基于Gauss解群优选的短弧初轨确定方法 | |
| CN114581532B (zh) | 多相机外参的联合标定方法、装置、设备和介质 | |
| US9709404B2 (en) | Iterative Kalman Smoother for robust 3D localization for vision-aided inertial navigation | |
| EP3842746A2 (en) | Dead reckoning method and apparatus for vehicle, device and storage medium | |
| WO2023004956A1 (zh) | 一种在高动态环境下的激光slam方法、系统设备及存储介质 | |
| EP3932781A1 (en) | Reverse trajectory tracking method and apparatus, electronic device and storage medium | |
| CN114013449A (zh) | 针对自动驾驶车辆的数据处理方法、装置和自动驾驶车辆 | |
| CN109612438B (zh) | 一种虚拟共面条件约束下的空间目标初轨确定方法 | |
| CN106909161B (zh) | 一种敏捷卫星零偏流角成像的最优姿态机动规划方法 | |
| US20150120096A1 (en) | Angles-Only Initial Orbit Determination (IOD) | |
| CN108692729A (zh) | 一种空间非合作目标相对导航协方差自适应修正滤波方法 | |
| CN108981750B (zh) | X射线脉冲双星光子序列仿真方法 | |
| CN108519110A (zh) | 基于图像信息的空间非合作目标自主相对导航在轨验证系统 | |
| CN112985391A (zh) | 一种基于惯性和双目视觉的多无人机协同导航方法和装置 | |
| CN111044037A (zh) | 一种光学卫星影像的几何定位方法及装置 | |
| CN115447584B (zh) | 一种车道中心线的确定方法、装置、设备及存储介质 | |
| CN108919283A (zh) | 一种星上自主的非合作目标相对导航方法和系统 | |
| CN105511483B (zh) | 鸟巢式星座及其设计方法 | |
| CN112070891A (zh) | 数字地面模型作为三维控制的影像区域网平差方法及系统 | |
| CN117029802A (zh) | 一种基于深度学习的多模态slam方法 | |
| CN120685066A (zh) | 一种基于李群的惯导/雷达/卫导三组合的发射系组合导航方法 | |
| WO2025242011A1 (zh) | 一种基于多源输入的三维模型重建方法和装置 | |
| JP6385380B2 (ja) | 演算装置、制御装置およびプログラム | |
| CN117346782B (zh) | 定位优化方法、装置、电子设备和存储介质 | |
| CN108759818A (zh) | 一种超高精度导星敏感器姿态确定的方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| 121 | Ep: the epo has been informed by wipo that ep was designated in this application |
Ref document number: 19919874 Country of ref document: EP Kind code of ref document: A1 |
|
| NENP | Non-entry into the national phase |
Ref country code: DE |
|
| 122 | Ep: pct application non-entry in european phase |
Ref document number: 19919874 Country of ref document: EP Kind code of ref document: A1 |
















