Patent Publication Number: US-11662209-B2

Title: Short arc initial orbit determining method based on gauss solution cluster

Description:
CROSS-REFERENCE TO RELATED APPLICATION 
     This application is a 371 of international application of PCT application serial no. PCT/CN2019/107186, filed on Sep. 23, 2019, which claims the priority benefit of China application no. 201910207587.6, filed on Mar. 19, 2019. The entirety of each of the above mentioned patent applications is hereby incorporated by reference herein and made a part of this specification. 
     BACKGROUND 
     Technical Field 
     The disclosure belongs to the field of astrodynamics, and more particularly relates to a preferable short arc initial orbit determining method based on Gauss solution cluster. 
     Description of Related Art 
     With the development of astrodynamics, the initial orbit determining technology has become one of the most important technologies nowadays. Conventional methods for determining the initial orbit are to calculate the orbits of satellites or space debris through Earth observation data. These methods normally adopt the equation of motions for solving problems between two objects. Due to the short observation time and less available data of the space-based observation system, the conventional method of determining the initial orbit is not suitable for space-based observation data. In addition, if the curvature of the track to be detected is too small, the initial orbit may not be determined. 
     The current methods for determining the initial orbit include: conventional Laplace method, Gauss method and related improved algorithms. The drawbacks of the conventional Gauss method are that too many roots are obtained after solving higher-order polynomials and need to be filtered; the improved Gauss method can avoid the problem of double roots, but it is necessary to choose a suitable initial value, otherwise the algorithm is calculated in multiple times and diverged or eventually a trivial solution or a wrong solution is obtained. 
     Current research data shows that there are many methods with excellent performance in the initial orbit determining field based on Earth observation data, and in the initial orbit determining field based on space-based observation data, the requirement for the time intervals of algorithms in related research is relatively long, and therefore the applicability of such short arc initial orbit determining methods is not high. The method of the disclosure can solve the problem of multiple root existing in the conventional Gauss method, and is applicable for the short arc initial orbit determining task. 
     SUMMARY 
     To overcome the shortcomings of the related art, the purpose of the present disclosure is to solve the technical problems of current technology, which has large errors in short arc initial orbit determining method based on conventional Gauss method and it is difficult to effectively choose a solution from multiple sets of solutions. 
     To achieve the above purpose, first of all, an embodiment of the present disclosure provides a preferable short arc initial orbit determining method based on a Gauss solution cluster, and the method includes the following steps: 
     S 1 . The observation data is divided into groups, and for each group of data, the Gauss method is adopted to find the target state vector for the corresponding time point to form a solution set of preliminary estimation. 
     S 2 . The solution set of preliminary estimation is divided into a position component vector solution set and a velocity component vector solution set, which are respectively clustered to obtain a position component vector solution cluster and a velocity component vector solution cluster. 
     S 3 . A two-dimensional trajectory solution set is generated based on the position component vector solution cluster and the velocity component vector solution cluster. 
     S 4 . The trajectory optimization method is adopted to evaluate each two-dimensional trajectory, calculate the number of roots of orbit corresponding to the optimal two-dimensional trajectory, and complete the determination of the initial orbit. 
     Specifically, step S 1  includes the following sub-steps: 
     S 101 . The observation data is divided into groups, and the vector data corresponding to every 3 time points is sorted into a group. 
     S 102 . For each group of vector data, the conventional Gauss method is adopted to find the target state vector at the corresponding time point to form a solution set of preliminary estimation. 
     Specifically, it is assumed that the observation time point is [t 0 ,t k ], and k+1≥3, and for input data [t 0 ,t k ] of a total of k+1 time points, the interval parameter Nrate=└rate*k┘ is taken, wherein └•┘ represents the rounding down, rate is the interval rate, and the grouped observation data is expressed as {{right arrow over (L)} t     s1   , {right arrow over (L)} t     s2   , {right arrow over (L)} t     s3   , {right arrow over (R)} t     s1   , {right arrow over (R)} t     s2   , {right arrow over (R)} t     s3   }, wherein {right arrow over (R)} t     s    represents a satellite observation position vector, {right arrow over (L)} t     s    represents an observation angle vector, 
     
       
         
           
             
               
                 ⌊ 
                 
                   
                     k 
                     - 
                     Nrate 
                   
                   2 
                 
                 ⌋ 
               
               ≤ 
               
                 s 
                 ⁢ 
                 2 
               
               ≤ 
               
                 ⌊ 
                 
                   
                     k 
                     - 
                     Nrate 
                   
                   2 
                 
                 ⌋ 
               
             
             , 
             
 
             
               
                 s 
                 ⁢ 
                 1 
               
               = 
               
                 ⌊ 
                 
                   
                     
                       2 
                       * 
                       s 
                       ⁢ 
                       2 
                     
                     + 
                     Nrate 
                     - 
                     k 
                     + 
                     1 
                   
                   2 
                 
                 ⌋ 
               
             
             , 
             
 
             
               
                 s 
                 ⁢ 
                 3 
               
               = 
               
                 
                   ⌊ 
                   
                     
                       
                         2 
                         * 
                         s 
                         ⁢ 
                         2 
                       
                       + 
                       Nrate 
                       - 
                       k 
                       + 
                       1 
                     
                     2 
                   
                   ⌋ 
                 
                 . 
               
             
           
         
       
     
     Specifically, step S 2  includes the following sub-steps: 
     S 201 . The non-reasonable solutions in the solution set of preliminary estimation are eliminated. 
     S 202 . The remaining state vector solution set is divided into a position component vector solution set and a velocity component vector solution set, and clustering is performed respectively so that correct solutions are gathered in the same cluster as many as possible, thereby obtaining the position component vector solution cluster and the velocity component vector solution cluster. 
     S 203 . Noise reduction preprocessing is performed on the position component vector solution cluster and the velocity component vector solution cluster, respectively. 
     Specifically, position data and velocity data satisfying any one of conditions |{right arrow over (r)} t     s   |&lt;6378.14 km or |{right arrow over (r)} t     s   &gt;42378 km or |{right arrow over (r)} t     s   &gt;11.2 km/s are excluded, wherein |{right arrow over (r)} t     s    represents the target position vector at time point t s , and |{right arrow over (v)} t     s    represents the target velocity vector at the time point t s . 
     Specifically, in step S 202 , the k-means clustering method is used for clustering and the Chauvenet&#39;s—criterion discriminating method is adopted to eliminate abnormal data. 
     Specifically, step S 3  includes the following sub-steps: 
     S 301 . The position component vector solution cluster and the velocity component vector solution cluster after the noise reduction are fitted to construct a state vector combination at various time points. 
     S 302 . A three-dimensional trajectory solution set of the target orbit is generated according to the state vector combination corresponding to various time points. 
     S 303 . The 3D trajectory solution set of the target orbit is projected to the instantaneous observation image plane according to the measurement status of the observation platform to obtain a 2D trajectory solution set. 
     Specifically, step S 4  includes the following sub-steps: 
     S 401 . The derivative error and position error for each two-dimensional trajectory are calculated. 
     The formula for calculating the derivative error of the m-th two-dimensional trajectory is as follows: 
     
       
         
           
             
               
                 Δ 
                 ⁢ 
                 
                   uv 
                   
                     m 
                     ⁢ 
                     ′ 
                   
                 
               
               = 
               
                 
                   ∑ 
                   
                     n 
                     = 
                     0 
                   
                   k 
                 
                   
                 
                   [ 
                   
                     
                       
                         ❘ 
                         &#34;\[LeftBracketingBar]&#34; 
                       
                       
                         
                           ( 
                           
                             
                               2 
                               ⁢ 
                               
                                 t 
                                 n 
                               
                               × 
                               
                                 
                                   U 
                                   ^ 
                                 
                                 
                                   ″ 
                                   ⁢ 
                                   m 
                                 
                               
                             
                             + 
                             
                               
                                 U 
                                 ^ 
                               
                               
                                 ′ 
                                 ⁢ 
                                 m 
                               
                             
                           
                           ) 
                         
                         - 
                         
                           ( 
                           
                             
                               2 
                               ⁢ 
                               
                                 t 
                                 n 
                               
                               × 
                               
                                 U 
                                 
                                   ″ 
                                   * 
                                 
                               
                             
                             + 
                             
                               U 
                               
                                 ′ 
                                 * 
                               
                             
                           
                           ) 
                         
                       
                       
                         ❘ 
                         &#34;\[RightBracketingBar]&#34; 
                       
                     
                     + 
                     
                       
                         ❘ 
                         &#34;\[LeftBracketingBar]&#34; 
                       
                       
                         
                           ( 
                           
                             
                               2 
                               ⁢ 
                               
                                 t 
                                 n 
                               
                               × 
                               
                                 
                                   V 
                                   ^ 
                                 
                                 
                                   ″ 
                                   ⁢ 
                                   m 
                                 
                               
                             
                             + 
                             
                               
                                 V 
                                 ^ 
                               
                               
                                 ′ 
                                 ⁢ 
                                 m 
                               
                             
                           
                           ) 
                         
                         - 
                         
                           ( 
                           
                             
                               2 
                               ⁢ 
                               
                                 t 
                                 n 
                               
                               × 
                               
                                 V 
                                 
                                   ″ 
                                   * 
                                 
                               
                             
                             + 
                             
                               V 
                               
                                 ′ 
                                 * 
                               
                             
                           
                           ) 
                         
                       
                       
                         ❘ 
                         &#34;\[RightBracketingBar]&#34; 
                       
                     
                   
                   ] 
                 
               
             
             ⁢ 
             
 
             
                                
               
                 
                   
                     U 
                     ^ 
                   
                   
                     t 
                     n 
                   
                   m 
                 
                 = 
                 
                   
                     
                       
                         U 
                         ^ 
                       
                       
                         ″ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       
                         U 
                         ^ 
                       
                       
                         ′ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     
                       U 
                       ^ 
                     
                     0 
                     m 
                   
                 
               
             
             ⁢ 
             
 
             
                                
               
                 
                   
                     V 
                     ^ 
                   
                   
                     t 
                     n 
                   
                   m 
                 
                 = 
                 
                   
                     
                       
                         V 
                         ^ 
                       
                       
                         ″ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       
                         V 
                         ^ 
                       
                       
                         ′ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     
                       V 
                       ^ 
                     
                     0 
                     m 
                   
                 
               
             
             ⁢ 
             
 
             
                               
               
                 
                   U 
                   
                     t 
                     n 
                   
                   * 
                 
                 = 
                 
                   
                     
                       U 
                       
                         ″ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       U 
                       
                         ′ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     U 
                     0 
                     * 
                   
                 
               
             
             ⁢ 
             
 
             
                               
               
                 
                   V 
                   
                     t 
                     n 
                   
                   * 
                 
                 = 
                 
                   
                     
                       V 
                       
                         ″ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       V 
                       
                         ′ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     V 
                     0 
                     * 
                   
                 
               
             
           
         
       
     
     wherein U″*, V″*, U′*, V′* are the quadratic polynomial fitting parameters of the two-dimensional trajectory {l*=(U t     n   , V t     n   *, t n , n=0, 1, 2, 3 . . . k)}, Û t     n     m {circumflex over (V)} t     n     m  is the corresponding image plane coordinates of the estimated trajectory l m ′ at the time point t n , m represents the m-th two-dimensional trajectory, and k is the total number of observation time points. 
     The formula for calculating the position error of the m-th two-dimensional trajectory is as follows: 
     
       
         
           
             
               
                 Rt 
                 m 
               
               = 
               
                 
                   ∑ 
                   
                     n 
                     = 
                     0 
                   
                   k 
                 
                   
                 
                   
                     D 
                     
                       t 
                       n 
                     
                     m 
                   
                   
                     Dc 
                     
                       t 
                       n 
                     
                     m 
                   
                 
               
             
             ⁢ 
             
 
             
               
                 D 
                 
                   t 
                   n 
                 
                 m 
               
               = 
               
                 
                   
                     
                       ( 
                       
                         
                           
                             U 
                             ^ 
                           
                           
                             t 
                             n 
                           
                           m 
                         
                         - 
                         
                           U 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                   + 
                   
                     
                       ( 
                       
                         
                           
                             V 
                             ^ 
                           
                           
                             t 
                             n 
                           
                           m 
                         
                         - 
                         
                           V 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                 
               
             
             ⁢ 
             
 
             
               
                 Dc 
                 
                   t 
                   n 
                 
                 m 
               
               = 
               
                 
                   
                     
                       ( 
                       
                         
                           U 
                           c 
                           m 
                         
                         - 
                         
                           U 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                   + 
                   
                     
                       ( 
                       
                         
                           V 
                           c 
                           m 
                         
                         - 
                         
                           V 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                 
               
             
           
         
       
     
     wherein, D t     s     m  is the image plane distance between the predicted value and the observed value, Dc t     s     m  is the image plane distance between the intersection point of the estimated trajectory and the observed trajectory and the observed point, to is the observation time, Û t     n     m {circumflex over (V)} t     n     m  is the corresponding image plane coordinate of the estimated trajectory l m ′ at the time point t n , U t     n   *V t     n   * is the corresponding image plane coordinate of the observed trajectory l* at the time point t n , and U c   m V c   m  is the image plane coordinate of the intersection point of the estimated trajectory l m ′ and the observed trajectory. 
     S 402 . The derivative error and the position error are taken into comprehensive consideration to select the optimal two-dimensional trajectory, and the number of root of orbit corresponding thereto is used as the optimal estimated value. 
     S 403 . The number of root of orbit corresponding to the optimal two-dimensional trajectory is output, and the determination of the initial orbit is completed. 
     Specifically, step S 402  is specifically as follows: 
     For the m-th orbit, the derivative error and position error are sorted from small to large, respectively, and the ranking sequence numbers Rank vel,m  and Rank dis,m  are obtained, the two sequence numbers are between 1˜N sv , and the sum of the two is calculated as the final error evaluation.
 
Rank m =Rank vel,m +Rank dis,m  
 
     For the final error evaluation, the smallest one m opt  is taken as the optimal predicted trajectory. 
     
       
         
           
             
               m 
               opt 
             
             = 
             
               
                 
                   
                     arg 
                     ⁢ 
                        
                     min 
                   
                      
                 
                 m 
               
               ⁢ 
               
                 { 
                 
                   
                     
                       Rank 
                       m 
                     
                     ; 
                     
                       m 
                       = 
                       1 
                     
                   
                   , 
                   2 
                   , 
                   
                     3 
                     ⁢ 
                         
                     … 
                     ⁢ 
                         
                     
                       N 
                       sv 
                     
                   
                 
                 } 
               
             
           
         
       
     
     wherein, N sv  is the number of orbit root obtained by the combination of the target position and the velocity data corresponding to various time points. 
     Second of all, an embodiment of the present disclosure provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the short arc initial orbit determining method described in the first aspect is implemented. 
     In general, compared with the related art, the above technical solution provided by the present disclosure has the following advantageous effects: 
     1. The present disclosure divides the observation data into groups and performs the conventional Gauss method on each group of data to increase the degree of data utilization. The number of root of short arc target orbits finally determined based on combination of multiple sets of data will be closer to the true value. 
     2. The present disclosure addresses the multiple root problem by removing all non-reasonable data and performing noise smoothing on all roots, calculating their corresponding orbit root, and adopting the trajectory-first method to conduct evaluation, finally the orbit root having the highest evaluation index is selected as the estimated solution. The method of the disclosure can effectively solve the multiple root problem existing in conventional Gauss method. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG.  1    is a flowchart of a short arc initial orbit determining method based on a Gauss solution cluster according to an embodiment of the present disclosure. 
         FIG.  2    is a schematic diagram showing space-based observation according to an embodiment of the present disclosure. 
         FIG.  3    is a solution cluster diagram of a target position obtained by the Gauss method according to an embodiment of the present disclosure. 
         FIG.  4    is a clustering diagram of k-means clustering performed on position solution cluster according to an embodiment of the present disclosure. 
         FIG.  5    is a position solution clustering diagram after the adopted Savitzky-Golay is smoothed according to an embodiment of the present disclosure. 
         FIG.  6    is a three-dimensional diagram corresponding to various position solutions obtained by calculating orbit roots according to an embodiment of the present disclosure. 
         FIG.  7    is a two-dimensional diagram showing projection of the three-dimensional orbit corresponding to various position solutions projected on the image plane according to an embodiment of the present disclosure. 
         FIG.  8    is a three-dimensional diagram of an optimal orbit and an actual observation orbit according to an embodiment of the present disclosure. 
         FIG.  9    is a flowchart of sub-steps of step S 2  based on a Gauss solution cluster according to an embodiment of the present disclosure. 
         FIG.  10    is a flowchart of sub-steps of step S 3  based on a Gauss solution cluster according to an embodiment of the present disclosure. 
         FIG.  11    is a flowchart of sub-steps of step S 4  based on a Gauss solution cluster according to an embodiment of the present disclosure. 
     
    
    
     DESCRIPTION OF THE EMBODIMENTS 
     In order to make the objectives, technical solutions, and advantages of the present disclosure clearer, the following further describes the present disclosure in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present disclosure and are not intended to limit the present disclosure. 
     First of all, some terms related to the present disclosure are explained as follows. 
     Solution Cluster: A Collection of Solutions 
     As shown in  FIG.  1   , the present disclosure provides a preferred short arc initial orbit determining method based on a Gauss solution cluster. The method includes the following steps: 
     S 1 . The observation data is divided into groups, and for each group of data, the Gauss method is adopted to find the target state vector for the corresponding time point to form a solution set of preliminary estimation. 
     S 2 . The solution set of preliminary estimation is divided into a position component vector solution set and a velocity component vector solution set, which are respectively clustered to obtain a position component vector solution cluster and a velocity component vector solution cluster. 
     S 3 . A two-dimensional trajectory solution set is generated based on the position component vector solution cluster and the velocity component vector solution cluster. 
     S 4 . The trajectory optimization method is adopted to evaluate each two-dimensional trajectory, calculate the number of root of orbit corresponding to the optimal two-dimensional trajectory, and complete the determination of the initial orbit. 
     Step S 1 . The observation data is divided into groups, and for each group of data, the Gauss method is adopted to find the target state vector for the corresponding time point to form a solution set of preliminary estimation. 
     S 101 . The observation data is divided into groups, and the vector data corresponding to every 3 time points is sorted into a group. 
     The observation data is observation data from space-based monocular imaging platform or Earth observation data. The observation data includes: the observation position (the position of the platform itself relative to the geocentric coordinate system) corresponding to each observation time point, the observation angle vector, satellite observation position vector represented by variant {right arrow over (R)} t     s    at time point t s , and observation angle vector represented by {right arrow over (T)} t     s    at time point t s . The position of the observation satellite, the angle of the observation camera, and the observation imaging data at the same time point constitute a pair of data, and multiple time point data constitute a group of data. 
     In order to construct a solution cluster subsequently, the present disclosure first needs to group the observation data, and the vector data corresponding to every 3 time points is sorted into a group. Assume that the observation time point is [t 0 ,t k ], and k+1≥3, and for input data [t 0 ,t k ] of a total of k+1 time points, the interval parameter Nrate=└rate*k┘ is taken, wherein └•┘ represents the rounding down, rate is the interval rate, normally the taken value is 0.2, and the grouped observation data is expressed as {{right arrow over (L)} t     s1   , {right arrow over (L)} t     s2   , {right arrow over (L)} t     s3   , {right arrow over (R)} t     s1   , {right arrow over (R)} t     s2   , {right arrow over (R)} t     s3   }, wherein {right arrow over (R)} t     s    represents a satellite observation position vector, {right arrow over (L)} t     s    represents an observation angle vector, 
     
       
         
           
             
               
                 ⌊ 
                 
                   
                     k 
                     - 
                     Nrate 
                   
                   2 
                 
                 ⌋ 
               
               ≤ 
               
                 s 
                 ⁢ 
                 2 
               
               ≤ 
               
                 ⌊ 
                 
                   
                     k 
                     - 
                     Nrate 
                   
                   2 
                 
                 ⌋ 
               
             
             , 
             
 
             
               
                 s 
                 ⁢ 
                 1 
               
               = 
               
                 ⌊ 
                 
                   
                     
                       2 
                       * 
                       s 
                       ⁢ 
                       2 
                     
                     + 
                     Nrate 
                     - 
                     k 
                     + 
                     1 
                   
                   2 
                 
                 ⌋ 
               
             
             , 
             
 
             
               
                 s 
                 ⁢ 
                 3 
               
               = 
               
                 
                   ⌊ 
                   
                     
                       
                         2 
                         * 
                         s 
                         ⁢ 
                         2 
                       
                       + 
                       Nrate 
                       - 
                       k 
                       + 
                       1 
                     
                     2 
                   
                   ⌋ 
                 
                 . 
               
             
           
         
       
     
     S 102 . For each group of vector data, the conventional Gauss method is adopted to find the target state vector at the corresponding time point to form a solution set of preliminary estimation. 
     Since the conventional Gauss method can only input three data and the input information is limited, the present disclosure pondered on using the Gauss method for multiple times to increase data utilization. The number of root of short arc target orbit finally determined based on a combination of multiple sets of data will be closer to the true value. 
     The Gauss method solves the target position vector and the target velocity vector at the intermediate time point s 2  by inputting the observation position vector and the observation angle vector corresponding to three different time points. The target refers to an unpowered target, which can be a debris or a runaway satellite, etc. According to the method of the present disclosure, every three time points in a series of time points are sorted into a group, and each group of data is solved by the conventional Gauss method once to obtain the target position vector and velocity vector solution at the intermediate time point in the group of data. The number of solutions obtained from each group of vector data may be 1, 2, and 3. The Gauss method is adopted multiple times to find the target vector corresponding to different time points to form a collection of solutions of preliminary estimation. 
     Step S 2 . The solution set of preliminary estimation is divided into a position component vector solution set and a velocity component vector solution set, which are respectively clustered to obtain a position component vector solution cluster and a velocity component vector solution cluster. 
     S 201 . The non-reasonable solutions in the solution set of preliminary estimation are eliminated. 
     The position data and velocity data satisfying any one of conditions |{right arrow over (r)} t     s   |&lt;6378.14 km (radius of Earth) or |{right arrow over (r)} t     s   |&gt;42378 km (geostationary orbit radius) or |{right arrow over (v)} t     s   |&gt;11.2 km/s (second cosmic velocity) are excluded, wherein {right arrow over (r)} t     s    represents the target position vector at time point t s , and {right arrow over (v)} t     s    represents the target velocity vector at the time point t s . 
     Obviously, a solution that satisfies any of the above three conditions is a non-reasonable solution. Therefore, the target position data in which the distance between the to-be-solved target position and center of Earth is smaller than the radius of the Earth or larger than the geostationary orbit radius and the target velocity data greater than the second cosmic velocity are excluded. 
     S 202 . The remaining state vector solution set is divided into a position component vector solution set and a velocity component vector solution set, and clustering is performed respectively so that correct solutions are gathered in the same cluster as many as possible, thereby obtaining the position component vector solution cluster and the velocity component vector solution cluster. 
     In order to separate the correct root from the wrong root so that the correct roots are gathered in the same cluster as many as possible, the position component vector solution set and the velocity component vector solution set are clustered respectively. The k-means clustering method is preferably adopted by the present disclosure for clustering (k is a value of 2 or 3), and Chauvenet&#39;s—criterion is adopted for discrimination to exclude abnormal data. The abnormal data refers to data that exhibits abnormality in the same group, such as outliers. 
     If the number of position solutions corresponding to various time points is less than or equal to 1, k-means clustering is not required; if the maximum number of position solutions corresponding to any time point is 2, then k-means clustering where k=2 is performed once; if the maximum number of position solutions corresponding to any time point is 3, a k-means clustering where k=3 is performed once. The velocity component vector solution set is clustered in the same way as described above. 
     S 203 . Noise reduction preprocessing is performed on the position component vector solution cluster and the velocity component vector solution cluster, respectively. 
     The present disclosure performs noise reduction preferably based on a Savitzky-Golay filtering method. 
     After all roots are processed by removing non-reasonable data and smoothing noise reduction, the accuracy of trajectory optimization can be improved. 
     Step S 3 . A two-dimensional trajectory solution set is generated based on the position component vector solution cluster and the velocity component vector solution cluster. 
     S 301 . The position component vector solution cluster and the velocity component vector solution cluster after the noise reduction are fitted to construct a state vector combination at various time points. 
     The data after noise deduction is utilized to perform fitting on the target position component vector and the target velocity component vector. The combination of velocity vector and position vector corresponding to the time point t s  is ({right arrow over (r)} t     s   ′, {right arrow over (v)} t     s   ′, t s ). 
     S 302 . A three-dimensional trajectory solution set of the target orbit is generated according to the state vector combination corresponding to various time points. 
     With the combination ({right arrow over (r)} t     s   ′, {right arrow over (v)} t     s   ′, t s ) of the velocity vector and position vector corresponding to a certain time point t s , the number of root {Ô m =(a m , e m , i m , Ω m , ω m , M m )} of orbit corresponding to ({right arrow over (r)} t     s   ′, {right arrow over (v)} t     s   ′, t s ) is calculated. Based on the number of root of orbit, it can be determined that the three-dimensional trajectory solution set is {{circumflex over (l)} m =({circumflex over (x)} t     n     m ,ŷ t     n     m ,{circumflex over (z)} t     n     m , t n )}, wherein a is the semi-major axis, e is the first eccentricity, i is the inclination of the orbit, Q is the right ascension of the ascending intersection, ω is the perigee angle, and M is the mean anomaly, ({circumflex over (x)} t     n     m ,{circumflex over (x)} t     n     m ,{circumflex over (x)} t     n     m ) is the corresponding three-dimensional coordinate of the observed trajectory {circumflex over (l)} m  at the time point t n , m=1, 2, . . . , N sv , n=0, 1, . . . , k, N sv  is the total number of root of orbit obtained based on the combination of the target positions and the velocity data corresponding to various time points. 
     S 303 . The 3D trajectory solution set of the target orbit is projected to the instantaneous observation image plane according to the measurement status of the observation platform to obtain a 2D trajectory solution set. 
     The 3D trajectory solution set {{circumflex over (l)} m =({circumflex over (x)} t     n     m ,ŷ t     n     m ,{circumflex over (z)} t     n     m ,t n )} of the target orbit is projected to the instantaneous observation image plane, thereby obtaining the 2D trajectory solution set {l m ′=Û t     n     m ,{circumflex over (V)} t     n     m , t n }, wherein {circumflex over (l)} m  is m-th 3D observed trajectories, ({circumflex over (x)} t     n     m ,{circumflex over (x)} t     n     m ,{circumflex over (x)} t     n     m ) is the corresponding 3D coordinate of the observed trajectory {circumflex over (l)} m  at the time point t n , l m ′ is the m-th estimated trajectory, Û t     n     m ,{circumflex over (V)} t     n     m  is the corresponding image plane coordinate of the estimated trajectory l m ′ at the time point t n , n=0, 1, . . . , k, m=1, 2, . . . , N sv . 
     Step S 4 . The trajectory optimization method is adopted to evaluate each two-dimensional trajectory, calculate the number of roots of orbit corresponding to the optimal two-dimensional trajectory, and complete the determination of the initial orbit. 
     S 401 . The derivative error and position error for each two-dimensional trajectory are calculated. 
     The formula for calculating the derivative error of the m-th two-dimensional trajectory is as follows: 
     
       
         
           
             
               
                 Δ 
                 ⁢ 
                 
                   uv 
                   
                     m 
                     ⁢ 
                     ′ 
                   
                 
               
               = 
               
                 
                   ∑ 
                   
                     n 
                     = 
                     0 
                   
                   k 
                 
                   
                 
                   [ 
                   
                     
                       
                         ❘ 
                         &#34;\[LeftBracketingBar]&#34; 
                       
                       
                         
                           ( 
                           
                             
                               2 
                               ⁢ 
                               
                                 t 
                                 n 
                               
                               × 
                               
                                 
                                   U 
                                   ^ 
                                 
                                 
                                   ″ 
                                   ⁢ 
                                   m 
                                 
                               
                             
                             + 
                             
                               
                                 U 
                                 ^ 
                               
                               
                                 ′ 
                                 ⁢ 
                                 m 
                               
                             
                           
                           ) 
                         
                         - 
                         
                           ( 
                           
                             
                               2 
                               ⁢ 
                               
                                 t 
                                 n 
                               
                               × 
                               
                                 U 
                                 
                                   ″ 
                                   * 
                                 
                               
                             
                             + 
                             
                               U 
                               
                                 ′ 
                                 * 
                               
                             
                           
                           ) 
                         
                       
                       
                         ❘ 
                         &#34;\[RightBracketingBar]&#34; 
                       
                     
                     + 
                     
                       
                         ❘ 
                         &#34;\[LeftBracketingBar]&#34; 
                       
                       
                         
                           ( 
                           
                             
                               2 
                               ⁢ 
                               
                                 t 
                                 n 
                               
                               × 
                               
                                 
                                   V 
                                   ^ 
                                 
                                 
                                   ″ 
                                   ⁢ 
                                   m 
                                 
                               
                             
                             + 
                             
                               
                                 V 
                                 ^ 
                               
                               
                                 ′ 
                                 ⁢ 
                                 m 
                               
                             
                           
                           ) 
                         
                         - 
                         
                           ( 
                           
                             
                               2 
                               ⁢ 
                               
                                 t 
                                 n 
                               
                               × 
                               
                                 V 
                                 
                                   ″ 
                                   * 
                                 
                               
                             
                             + 
                             
                               V 
                               
                                 ′ 
                                 * 
                               
                             
                           
                           ) 
                         
                       
                       
                         ❘ 
                         &#34;\[RightBracketingBar]&#34; 
                       
                     
                   
                   ] 
                 
               
             
             ⁢ 
             
 
             
                                
               
                 
                   
                     U 
                     ^ 
                   
                   
                     t 
                     n 
                   
                   m 
                 
                 = 
                 
                   
                     
                       
                         U 
                         ^ 
                       
                       
                         ″ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       
                         U 
                         ^ 
                       
                       
                         ′ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     
                       U 
                       ^ 
                     
                     0 
                     m 
                   
                 
               
             
             ⁢ 
             
 
             
                                
               
                 
                   
                     V 
                     ^ 
                   
                   
                     t 
                     n 
                   
                   m 
                 
                 = 
                 
                   
                     
                       
                         V 
                         ^ 
                       
                       
                         ″ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       
                         V 
                         ^ 
                       
                       
                         ′ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     
                       V 
                       ^ 
                     
                     0 
                     m 
                   
                 
               
             
             ⁢ 
             
 
             
                                
               
                 
                   U 
                   
                     t 
                     n 
                   
                   * 
                 
                 = 
                 
                   
                     
                       U 
                       
                         ″ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       U 
                       
                         ′ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     U 
                     0 
                     * 
                   
                 
               
             
             ⁢ 
             
 
             
                                
               
                 
                   V 
                   
                     t 
                     n 
                   
                   * 
                 
                 = 
                 
                   
                     
                       V 
                       
                         ″ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       V 
                       
                         ′ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     V 
                     0 
                     * 
                   
                 
               
             
           
         
       
     
     wherein U″*, V″*, U′*, V′* are the quadratic polynomial fitting parameters of the two-dimensional trajectory {l*=(U t     n   , V t     n   *, t n , n=0, 1, 2, 3 . . . k)}, Û t     n     m ,{circumflex over (V)} t     n     m  is the corresponding image plane coordinates of the estimated trajectory l m ′ at the time point t n , m represents the m-th two-dimensional trajectory, and k is the total number of observation time points. 
     The camera parameters (orientation) remain the same, and the phase plane coordinates (U t     n   *, V t     n   *) of various time points with respect to the observation of target can be obtained. Quadratic polynomial fitting is performed on the observation data for two times in U and V directions respectively, thereby obtaining the second term and the fourth term. 
     The formula for calculating the position error of the m-th two-dimensional trajectory is as follows: 
     
       
         
           
             
               
                 Rt 
                 m 
               
               = 
               
                 
                   ∑ 
                   
                     n 
                     = 
                     0 
                   
                   k 
                 
                   
                 
                   
                     D 
                     
                       t 
                       n 
                     
                     m 
                   
                   
                     Dc 
                     
                       t 
                       n 
                     
                     m 
                   
                 
               
             
             ⁢ 
             
 
             
               
                 D 
                 
                   t 
                   n 
                 
                 m 
               
               = 
               
                 
                   
                     
                       ( 
                       
                         
                           
                             U 
                             ^ 
                           
                           
                             t 
                             n 
                           
                           m 
                         
                         - 
                         
                           U 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                   + 
                   
                     
                       ( 
                       
                         
                           
                             V 
                             ^ 
                           
                           
                             t 
                             n 
                           
                           m 
                         
                         - 
                         
                           V 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                 
               
             
             ⁢ 
             
 
             
               
                 Dc 
                 
                   t 
                   n 
                 
                 m 
               
               = 
               
                 
                   
                     
                       ( 
                       
                         
                           U 
                           c 
                           m 
                         
                         - 
                         
                           U 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                   + 
                   
                     
                       ( 
                       
                         
                           V 
                           c 
                           m 
                         
                         - 
                         
                           V 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                 
               
             
           
         
       
     
     wherein, D t     n     m  is the image plane distance between the predicted value and the observed value, Dc t     n     m  is the image plane distance between the intersection point of the estimated trajectory and the observed trajectory and the observed point, t n  is the observation time, Û t     n     m ,{circumflex over (V)} t     n     m  is the corresponding image plane coordinate of the estimated trajectory l m ′ at the time point t n , U t     n   *V t     n   * is the corresponding image plane coordinate of the observed trajectory l* at the time point t n , and U c   m V c   m  is the image plane coordinate of the intersection point of the estimated trajectory l m ′ and the observed trajectory. 
     S 402 . The derivative error and the position error are taken into comprehensive consideration to select the optimal two-dimensional trajectory, and the number of root of orbit corresponding thereto is used as the optimal estimated value. 
     For the m-th orbit, the derivative error and position error are sorted from small to large, respectively, and the ranking sequence numbers Rank vel,m  and Rank dis,m  are obtained, the two sequence numbers are between 1˜N sv , and the sum of the two is calculated as the final error evaluation.
 
Rank m =Rank vel,m +Rank dis,m  
 
     For the final error evaluation, take the smallest value m opt  as the optimal two-dimensional trajectory: 
     
       
         
           
             
               m 
               opt 
             
             = 
             
               
                 
                   
                     arg 
                     ⁢ 
                        
                     min 
                   
                      
                 
                 m 
               
               ⁢ 
               
                 { 
                 
                   
                     
                       Rank 
                       m 
                     
                     ; 
                     
                       m 
                       = 
                       1 
                     
                   
                   , 
                   2 
                   , 
                   
                     3 
                     ⁢ 
                         
                     … 
                     ⁢ 
                         
                     
                       N 
                       sv 
                     
                   
                 
                 } 
               
             
           
         
       
     
     wherein, N sv  is the number of root of orbits obtained based on the combination of the target position and the velocity data corresponding to various time points. 
     S 403 . The number of root of orbit corresponding to the optimal two-dimensional trajectory is output, and the determination of the initial orbit is completed. 
     The present disclosure is suitable for single space-based imaging observation platform, long-range unpowered targets under short-term observation conditions (greater than 5 s), such as debris, out-of-control satellites, etc. The disclosure can effectively solve the problem of multiple roots that cannot be solved by the conventional Gauss method, and has higher utilization of original observation data. Moreover, the present disclosure adopts the concept of solution cluster and optimal method to improve the accuracy of the initial orbit determination. Under the same input error conditions, the technical solution of the present disclosure has better performance than other schemes in terms of the short-term observation data, that is, the present disclosure is more suitable for short-arc observation than other disclosures. 
     Example 1 
     As shown in Table 1, Table 1 shows data of observation satellites and number of root of target orbits in Example 1. 
     
       
         
           
               
               
               
               
               
               
             
               
                 TABLE 1 
               
               
                   
               
               
                   
                 a (km) 
                 e 
                 i (°) 
                 u (°) 
                 Ω (°) 
               
               
                   
               
             
            
               
                   
               
            
           
           
               
               
               
               
               
               
            
               
                 Observing 
                 6978.1 
                 0 
                 97.6893 
                 0.1719 
                 275.5927 
               
               
                 space- 
                   
                   
                   
                   
                   
               
               
                 based 
                   
                   
                   
                   
                   
               
               
                 data 
                   
                   
                   
                   
                   
               
               
                 Target 
                 6987.0 
                 7.8e−09 
                 79.1255 
                 4.7555 
                 272.7279 
               
               
                 data 
               
               
                   
               
            
           
         
       
     
     S 101 . The observation data is divided into groups, and the vector data corresponding to every 3 time points is sorted into a group. 
     With respect to the input data [t 0 , t 170 ] of a total of 171 time points, the interval parameter Nrate=└rate*k┘ is taken, the observation data {{right arrow over (L)} t     s1   , {right arrow over (L)} t     s2   , {right arrow over (L)} t     s3   , {right arrow over (R)} t     s1   , {right arrow over (R)} t     s2   , {right arrow over (R)} t     s3   } is sorted into a group, wherein └•┘ represents the rounding down, rate is the interval rate, and 
               rate   =   0.2     ,       s   ⁢   1     =     ⌊         2   *   s   ⁢   2     +   34   -   170   +   1     2     ⌋       ,   
       ⌊       170   -   34     2     ⌋     ≤     s   ⁢   2     ≤     ⌊       170   +   34     2     ⌋       ,   
       s   ⁢   3     =     ⌊         2   *   s   ⁢   2     -   34   +   170   -   1     2     ⌋       ,         
{right arrow over (R)} t     s    represents a satellite observation position vector, {right arrow over (L)} t     s    represents an observation angle vector.
 
     With respect to the observation satellites and the number of root of target orbits in the data in Table 1, simulation is conducted to obtain the data regarding observation satellites and partial position of the target as shown in Table 2. 
     
       
         
           
               
               
               
               
               
             
               
                 TABLE 2 
               
             
            
               
                   
               
               
                 Point 
                   
                 {right arrow over (r)} t  (km) 
                 {circumflex over (L)} t  (incl. error/km) 
                 {right arrow over (R)} t  (km) 
               
            
           
           
               
               
               
               
               
               
               
               
               
               
               
            
               
                 number 
                 t (s) 
                 x t     s     
                 y t     s     
                 z t     s     
                 x t     s     
                 y t     s     
                 z t     s     
                 x t     s     
                 y t     s     
                 z t     s     
               
               
                   
               
            
           
           
               
               
               
               
               
               
               
               
               
               
               
            
               
                 0 
                 0.0 
                 306.30 
                 −6980.14 
                 −52.03 
                 −0.9877 
                 −0.0934 
                 −0.1254 
                 705.00 
                 −6942.43 
                 −1.41 
               
            
           
           
               
               
               
               
               
            
               
                   
                   
                 . . . 
                 . . . 
                 . . . 
               
            
           
           
               
               
               
               
               
               
               
               
               
               
               
            
               
                 34 
                 3.4 
                 311.15 
                 −6980.06 
                 −26.81 
                 −0.9872 
                 −0.09439 
                 −0.1285 
                 701.58 
                 −6942.73 
                 24.01 
               
            
           
           
               
               
               
               
               
            
               
                   
                   
                 . . . 
                 . . . 
                 . . . 
               
            
           
           
               
               
               
               
               
               
               
               
               
               
               
            
               
                 68 
                 6.8 
                 315.99 
                 −6979.90 
                 −1.59 
                 −0.9867 
                 −0.0954 
                 −0.1318 
                 698.15 
                 −6942.94 
                 49.44 
               
            
           
           
               
               
               
               
               
            
               
                   
                   
                 . . . 
                 . . . 
                 . . . 
               
            
           
           
               
               
               
               
               
               
               
               
               
               
               
            
               
                 102 
                 10.2 
                 320.84 
                 −6979.64 
                 23.62 
                 −0.9861 
                 −0.0964 
                 −0.1351 
                 694.71 
                 −6943.06 
                 74.87 
               
            
           
           
               
               
               
               
               
            
               
                   
                   
                 . . . 
                 . . . 
                 . . . 
               
            
           
           
               
               
               
               
               
               
               
               
               
               
               
            
               
                 136 
                 13.6 
                 325.67 
                 −6979.28 
                 48.84 
                 −0.9855 
                 −0.0976 
                 −0.1387 
                 691.26 
                 −6943.09 
                 100.30 
               
            
           
           
               
               
               
               
               
            
               
                   
                   
                 . . . 
                 . . . 
                 . . . 
               
            
           
           
               
               
               
               
               
               
               
               
               
               
               
            
               
                 170 
                 17.0 
                 330.37 
                 −6978.85 
                 73.31 
                 −0.9848 
                 −0.0987 
                 −0.1424 
                 687.80 
                 −6943.01 
                 125.73 
               
               
                   
               
            
           
         
       
     
     S 102 . For each group of vector data, the conventional Gauss method is adopted to find the target state vector at the corresponding time point to form a solution set of preliminary estimation. 
     A single group of data after grouping is as shown in  FIG.  2   , point E represents the center of the coordinate system. The distance between the observation point and the target is defined as ρ i , then: 
     
       
         
           
             { 
             
               
                 
                   
                     
                       
                         
                           r 
                           1 
                         
                         → 
                       
                       = 
                       
                         
                           
                             R 
                             1 
                           
                           → 
                         
                         + 
                         
                           
                             ρ 
                             1 
                           
                           ⁢ 
                           
                             
                               L 
                               1 
                             
                             → 
                           
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           r 
                           2 
                         
                         → 
                       
                       = 
                       
                         
                           
                             R 
                             2 
                           
                           → 
                         
                         + 
                         
                           
                             ρ 
                             2 
                           
                           ⁢ 
                           
                             
                               L 
                               2 
                             
                             → 
                           
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           r 
                           3 
                         
                         → 
                       
                       = 
                       
                         
                           
                             R 
                             3 
                           
                           → 
                         
                         + 
                         
                           
                             ρ 
                             3 
                           
                           ⁢ 
                           
                             
                               L 
                               3 
                             
                             → 
                           
                         
                       
                     
                   
                 
               
                
             
           
         
       
     
     wherein, the satellite observation position vector at the time point t i  and the corresponding position vector of the target are {right arrow over (R)} i  and {right arrow over (r)} i  respectively, and the unit sight vector corresponding to the observation time points t 1 , t 2  and t 3  are {right arrow over (L)} 1 , {right arrow over (L)} 2  and {right arrow over (L)} 3  respectively. There are six unknown vectors in the above formula, namely, {right arrow over (r)} 1 , {right arrow over (r)} 2  and {right arrow over (r)} 3  as well as ρ 1 , ρ 2  and ρ 3 ; r 2  is the magnitude of {right arrow over (r)} 2 , and the obtained target position is as shown in  FIG.  3   , the target velocity is as shown in  FIG.  4   , wherein the mark * indicates the calculated position solution, the mark + indicates the actual position solution of target, and the mark ▴ indicates the observation satellite position. According to the Gauss method, the initial value of {right arrow over (r)} 2  and {right arrow over (v)} 2  can be solved, and the data obtained from this example is shown in Table 3. Table 3 shows data concerning target position and velocity solution cluster based on Gauss method in Example 1. 
     
       
         
           
               
               
               
               
               
             
               
                 TABLE 3 
               
               
                   
               
               
                   
                   
                 position vector 
                 velocity vector 
                 Position vector data 
               
               
                 Combination 
                 Point t     s1     
                 solution set based 
                 solution set based on 
                 after removal of non- 
               
               
                 of 
                 Point t     s2     
                 on Gauss method 
                 Gauss method (m/s, 
                 reasonable data (km, 
               
               
                 sequence 
                 Point t     s3     
                 (km, km, km) 
                 m/s, m/s) 
                 km, km) 
               
               
                   
               
             
            
               
                  1 
                 Point 0, 
                 (−162380.97, 
                 (1037557.63, 47310.82, 
                 (256.55, −6985.65, 
               
               
                   
                 Point 68, 
                 −22712.67, 
                 −19460.97); 
                 −9.53); 
               
               
                   
                 Point 
                 −21732.39); 
                 (1801.66, 79.89,  
                 (880.40, −6925.32, 
               
               
                   
                 136 
                 (256.55, 
                 7407.1059); 
                 73.78) 
               
               
                   
                   
                 −6985.65, −9.53); 
                 (−2171.29, −101.28, 
                   
               
               
                   
                   
                 (880.40, 
                 7510.17) 
                   
               
               
                   
                   
                 −6925.32, 73.78) 
                   
                   
               
               
                 . . . 
                 . . . 
                 . . . 
                 . . . 
                 . . . 
               
               
                 34 
                 Point 
                 (−145427.83, 
                 (948881.44, 43250.12, 
                 (247.73, −6986.78, 
               
               
                   
                 34, 
                 −21233.37, 
                 −17276.82); 
                 12.88); 
               
               
                   
                 Point 
                 −19943.24); 
                 (1892.90, 110.84, 
                 (874.68, −6925.47, 
               
               
                   
                 102, 
                 (247.73, −6986.78, 
                 7404.08); 
                 98.77) 
               
               
                   
                 Point 
                 12.88); 
                 (−2182.56, 
                   
               
               
                   
                 170 
                 (874.68, −6925.47, 
                 −74.81, 7510.29) 
                   
               
               
                   
                   
                 98.77) 
               
               
                   
               
            
           
         
       
     
     S 201 . The non-reasonable solutions in the solution set of preliminary estimation are eliminated. 
     The principle for elimination is that positions and velocity data that satisfy |{right arrow over (r)} t     s   |&lt;6378.14 km, |{right arrow over (r)} t     s   |&gt;42378 km and |{right arrow over (v)} t     s   |&gt;11.2 km/s are excluded. The obtained positions and velocity in the example are as shown in Table 4, which shows the target position error data before smoothing in Example 1. 
     
       
         
           
               
               
               
               
               
             
               
                 TABLE 4 
               
               
                   
               
               
                   
                   
                 Optimized  
                 Simulated  
                 Position 
               
               
                   
                 Point t     s1     
                 position 
                 position 
                 vector  
               
               
                   
                 Point t     s2     
                 vector  
                 vector  
                 error 
               
               
                 Set 
                 Point t     s3     
                 (km, km, km) 
                 (km, km, km) 
                 (km, km, km) 
               
               
                   
               
             
            
               
                  1 
                 Point 0, 
                 (256.55,  
                 (315.99, 
                 (59.44, 5.74, 
               
               
                   
                 Point 68, 
                 −6985.65, 
                 −6979.90, 
                 7.94) 
               
               
                   
                 Point 136 
                 −9.53) 
                 −1.59) 
                 (−564.40, −54.57, 
               
               
                   
                   
                 (880.40, 
                   
                 −75.38) 
               
               
                   
                   
                 −6925.32, 
                   
                   
               
               
                   
                   
                 73.78) 
                   
                   
               
               
                 . . . 
                 . . . 
                 . . . 
                 . . . 
                 . . . 
               
               
                 34 
                 Point 34, 
                 (247.73, 
                 (320.84, 
                 (73.10, 7.14, 
               
               
                   
                 Point 102, 
                 −6986.78,  
                 −6979.64,  
                 10.73) 
               
               
                   
                 Point 170 
                 12.88) 
                 23.62) 
                 (−553.84, −54.16, 
               
               
                   
                   
                 (874.68, 
                   
                 −75.14) 
               
               
                   
                   
                 −6925.47,  
                   
                   
               
               
                   
                   
                 98.77) 
               
               
                   
               
            
           
         
       
     
     It can be expressed as follows. 
     
       
         
           
             
               { 
               
                 
                   
                     ( 
                     
                       
                         
                           r 
                           → 
                         
                         
                           t 
                           s 
                         
                         m 
                       
                       , 
                       
                         t 
                         s 
                       
                     
                     ) 
                   
                   ; 
                   
                     m 
                     = 
                     1 
                   
                 
                 , 
                 2 
                 , 
                 
                   
                     3 
                     ⁢ 
                         
                     … 
                     ⁢ 
                         
                     N 
                   
                   ; 
                   
                     ( 
                     
                       
                         ⌊ 
                         
                           
                             170 
                             - 
                             34 
                           
                           2 
                         
                         ⌋ 
                       
                       ≤ 
                       s 
                       ≤ 
                       
                         ⌊ 
                         
                           
                             170 
                             + 
                             34 
                           
                           2 
                         
                         ⌋ 
                       
                     
                     ) 
                   
                 
               
               } 
             
             ⁢ 
             
 
             
               { 
               
                 
                   
                     ( 
                     
                       
                         
                           v 
                           → 
                         
                         
                           t 
                           s 
                         
                         m 
                       
                       , 
                       
                         t 
                         s 
                       
                     
                     ) 
                   
                   ; 
                   
                     m 
                     = 
                     1 
                   
                 
                 , 
                 2 
                 , 
                 
                   
                     3 
                     ⁢ 
                         
                     … 
                     ⁢ 
                         
                     N 
                   
                   ; 
                   
                     ( 
                     
                       
                         ⌊ 
                         
                           
                             170 
                             - 
                             34 
                           
                           2 
                         
                         ⌋ 
                       
                       ≤ 
                       s 
                       ≤ 
                       
                         ⌊ 
                         
                           
                             170 
                             + 
                             34 
                           
                           2 
                         
                         ⌋ 
                       
                     
                     ) 
                   
                 
               
               } 
             
           
         
       
     
     wherein, N is the number of remaining positions or velocity solutions after removing non-reasonable data. 
     S 202 . The remaining state vector solution set is divided into a position component vector solution set and a velocity component vector solution set, and clustering is performed respectively so that correct solutions are gathered in the same cluster as many as possible, thereby obtaining the position component vector solution cluster and the velocity component vector solution cluster. 
     Perform k-means clustering on the obtained 
             {       (         r   →       t   s     m     ,     t   s       )     ,       (         v   →       t   s     m     ,     t   s       )     ;     m   =   1       ,   2   ,       3   ⁢         …   ⁢         N     ;     (       ⌊       170   -   34     2     ⌋     ≤   s   ≤     ⌊       170   +   34     2     ⌋       )         }         
respectively. Preferably, if the number of position solutions corresponding to all the time points is 1, k-means clustering is not required; if the maximum number of position solutions corresponding to any time point is 2, then k-means clustering where k=2 is performed once; if the maximum number of position solutions corresponding to any time point is 3, a k-means clustering where k=3 is performed once. The velocity component vector solution set is clustered in the same way as described above.
 
     Taking the location solution as an example, the clustering result is shown in  FIG.  5   , wherein the mark * indicates cluster  1 , the mark o indicates cluster  2 , the mark + is the target actual positions, and the mark ▴ indicates the observation satellite position. The data in this example is as shown in Table 5, which shows the target position error data in Example 1 after smoothing. 
     
       
         
           
               
               
               
               
               
               
               
               
             
               
                 TABLE 5 
               
               
                   
               
               
                   
                   a  (km) 
                 
                   e 
                 
                 incl (°) 
                   ω  (°) 
                 Ω (°) 
                 M (°) 
                   u  (°) 
               
               
                   
               
             
            
               
                   
               
            
           
           
               
               
               
               
               
               
               
               
            
               
                 Target 
                 6987.0 
                 7.8e−09 
                 79.1255 
                 / 
                 272.7279 
                 / 
                 4.7555 
               
               
                 data 
                   
                   
                   
                   
                   
                   
                   
               
               
                 Data of 
                 6928.366 
                 0.008 
                 80.6249 
                 185.7025 
                 272.7467 
                 173.7632 
                 −0.5343 
               
               
                 present 
                   
                   
                   
                   
                   
                   
                   
               
               
                 method 
                   
                   
                   
                   
                   
                   
                   
               
               
                 Original 
                 −0.236 
                 143393.4 
                 84.1812 
                 88.1617 
                 2.7930 
                 0 
                 88.1617 
               
               
                 solution 
                   
                   
                   
                   
                   
                   
                   
               
               
                 data 1 
                   
                   
                   
                   
                   
                   
                   
               
               
                 based on 
                   
                   
                   
                   
                   
                   
                   
               
               
                 Gauss 
                   
                   
                   
                   
                   
                   
                   
               
               
                 method 
                   
                   
                   
                   
                   
                   
                   
               
               
                 Original 
                 8447.614 
                 0.1748 
                 112.7304 
                 11.4240 
                 278.8715 
                 352.7453 
                 4.1693 
               
               
                 solution 
                   
                   
                   
                   
                   
                   
                   
               
               
                 data 2 
                   
                   
                   
                   
                   
                   
                   
               
               
                 based on 
                   
                   
                   
                   
                   
                   
                   
               
               
                 Gauss 
                   
                   
                   
                   
                   
                   
                   
               
               
                 method 
                   
                   
                   
                   
                   
                   
                   
               
               
                 Original 
                 6805.225 
                 0.0262 
                 86.0354 
                 189.0760 
                 273.7521 
                 170.6514 
                 −0.2726 
               
               
                 solution 
                   
                   
                   
                   
                   
                   
                   
               
               
                 data 3 
                   
                   
                   
                   
                   
                   
                   
               
               
                 based on 
                   
                   
                   
                   
                   
                   
                   
               
               
                 Gauss 
                   
                   
                   
                   
                   
                   
                   
               
               
                 method 
               
               
                   
               
            
           
         
       
     
     Example 1 adopts Chauvenet&#39;s —criterion for discrimination to eliminate abnormal data. 
     When the confidence probability is taken as 
               1   -     1     2   ⁢   n         ,         
the confidence probability formula of Chauvenet&#39;s—criterion can be expressed as: w n =1+0.4In(n).
 
     With respect to the value x i  to be tested, the absolute value of the difference between the value x i  and the sample mean is calculated. If the following formula is satisfied, the current value to be tested is excluded: |x i − x |&gt;w n S x , wherein x i  is the value to be tested. Here, the components in the three directions (x, y, z) of the position and velocity vectors are selected for three detections respectively. If any one of the detections satisfies the above formula, then it is excluded.  x  is the average value of such sample, w n  is the confidence probability of the Chauvenet&#39;s criterion, S x  is the standard deviation of such sample. 
     S 203 . Noise reduction preprocessing is performed on the position component vector solution cluster and the velocity component vector solution cluster, respectively. 
     The present disclosure performs noise reduction preferably based on a Savitzky-Golay filtering method, and the processing is as follows: 
     
       
         
           
             
               
                 ( 
                 
                   
                     
                       
                         
                           x 
                           _ 
                         
                         
                           t 
                           
                             s 
                             - 
                             2 
                           
                         
                       
                     
                     
                       
                         
                           x 
                           _ 
                         
                         
                           t 
                           
                             s 
                             - 
                             1 
                           
                         
                       
                     
                     
                       
                         
                           x 
                           _ 
                         
                         
                           t 
                           s 
                         
                       
                     
                     
                       
                         
                           x 
                           _ 
                         
                         
                           t 
                           
                             s 
                             + 
                             1 
                           
                         
                       
                     
                     
                       
                         
                           x 
                           _ 
                         
                         
                           t 
                           
                             s 
                             + 
                             1 
                           
                         
                       
                     
                   
                   
                     
                       
                         
                           y 
                           _ 
                         
                         
                           t 
                           
                             s 
                             - 
                             1 
                           
                         
                       
                     
                     
                       
                         
                           y 
                           _ 
                         
                         
                           t 
                           
                             s 
                             - 
                             1 
                           
                         
                       
                     
                     
                       
                         
                           y 
                           _ 
                         
                         
                           t 
                           s 
                         
                       
                     
                     
                       
                         
                           y 
                           _ 
                         
                         
                           t 
                           
                             s 
                             + 
                             1 
                           
                         
                       
                     
                     
                       
                         
                           y 
                           _ 
                         
                         
                           t 
                           
                             s 
                             + 
                             1 
                           
                         
                       
                     
                   
                   
                     
                       
                         
                           
                             
                               z 
                                 
                             
                             _ 
                           
                             
                         
                         
                           t 
                           
                             s 
                             - 
                             1 
                           
                         
                       
                     
                     
                       
                         
                           
                             z 
                             _ 
                           
                             
                         
                         
                           t 
                           
                             s 
                             - 
                             1 
                           
                         
                       
                     
                     
                       
                         
                           
                             z 
                             _ 
                           
                             
                         
                         
                           t 
                           s 
                         
                       
                     
                     
                       
                         
                           
                             z 
                             _ 
                           
                             
                         
                         
                           t 
                           
                             s 
                             + 
                             1 
                           
                         
                       
                     
                     
                       
                         
                           
                             z 
                             _ 
                           
                             
                         
                         
                           t 
                           
                             s 
                             + 
                             1 
                           
                         
                       
                     
                   
                 
                 ) 
               
               T 
             
             = 
             
               
                 1 
                 70 
               
               ⁢ 
               
                 ( 
                 
                   
                     
                       69 
                     
                     
                       4 
                     
                     
                       
                         - 
                         6 
                       
                     
                     
                       4 
                     
                     
                       
                         - 
                         1 
                       
                     
                   
                   
                     
                       4 
                     
                     
                       54 
                     
                     
                       24 
                     
                     
                       
                         - 
                         16 
                       
                     
                     
                       4 
                     
                   
                   
                     
                       
                         - 
                         6 
                       
                     
                     
                       24 
                     
                     
                       34 
                     
                     
                       24 
                     
                     
                       
                         - 
                         6 
                       
                     
                   
                   
                     
                       4 
                     
                     
                       
                         - 
                         16 
                       
                     
                     
                       24 
                     
                     
                       54 
                     
                     
                       4 
                     
                   
                   
                     
                       
                         - 
                         1 
                       
                     
                     
                       4 
                     
                     
                       
                         - 
                         6 
                       
                     
                     
                       4 
                     
                     
                       69 
                     
                   
                 
                 ) 
               
               ⁢ 
               
                 
                   ( 
                   
                     
                       
                         
                           x 
                           
                             t 
                             
                               s 
                               - 
                               2 
                             
                           
                         
                       
                       
                         
                           x 
                           
                             t 
                             
                               s 
                               - 
                               1 
                             
                           
                         
                       
                       
                         
                           x 
                           
                             t 
                             s 
                           
                         
                       
                       
                         
                           x 
                           
                             t 
                             
                               s 
                               + 
                               1 
                             
                           
                         
                       
                       
                         
                           x 
                           
                             t 
                             
                               s 
                               + 
                               2 
                             
                           
                         
                       
                     
                     
                       
                         
                           y 
                           
                             t 
                             
                               s 
                               - 
                               2 
                             
                           
                         
                       
                       
                         
                           y 
                           
                             t 
                             
                               s 
                               - 
                               1 
                             
                           
                         
                       
                       
                         
                           y 
                           
                             t 
                             s 
                           
                         
                       
                       
                         
                           y 
                           
                             t 
                             
                               s 
                               + 
                               1 
                             
                           
                         
                       
                       
                         
                           y 
                           
                             t 
                             
                               s 
                               + 
                               2 
                             
                           
                         
                       
                     
                     
                       
                         
                           z 
                           
                             t 
                             
                               s 
                               - 
                               2 
                             
                           
                         
                       
                       
                         
                           z 
                           
                             t 
                             
                               s 
                               - 
                               1 
                             
                           
                         
                       
                       
                         
                           z 
                           
                             t 
                             s 
                           
                         
                       
                       
                         
                           z 
                           
                             t 
                             
                               s 
                               + 
                               1 
                             
                           
                         
                       
                       
                         
                           z 
                           
                             t 
                             
                               s 
                               + 
                               2 
                             
                           
                         
                       
                     
                   
                   ) 
                 
                 T 
               
             
           
         
       
     
     After using the above-mentioned five-point and three-order template for smoothing, the noise can be effectively reduced in three directions (x, y, z). The position data after smoothing will be more concentrated in the solution space. Take the position solution cluster as an example, the results before and after smoothing are shown in  FIG.  6   , wherein the mark * indicates the target position before smoothing, the mark o indicates the target position after smoothing. 
     S 301 . The position component vector solution cluster and the velocity component vector solution cluster after the noise reduction are fitted to construct a state vector combination at various time points. 
     For the data after smoothing and noise-reduction, a quadratic polynomial fitting is performed in three directions of x, y, and z. After calculating the estimated values of the target positions corresponding to each intermediate time point 
             {       t   s   ′     ;       ⌊       170   -   34     2     ⌋     ≤   s   ≤     ⌊       170   +   34     2     ⌋         }         
on the x, y, and z components, they can be combined to obtain the corresponding position correction estimated value
 
             {           r   →       t   s       ′     ;       ⌊       170   -   34     2     ⌋     ≤   s   ≤     ⌊       170   +   34     2     ⌋         }         
and velocity correction estimated value
 
     
       
         
           
             
               { 
               
                 
                   
                     
                       r 
                       → 
                     
                     
                       t 
                       s 
                     
                   
                   ′ 
                 
                 ; 
                 
                   
                     ⌊ 
                     
                       
                         170 
                         - 
                         34 
                       
                       2 
                     
                     ⌋ 
                   
                   ≤ 
                   s 
                   ≤ 
                   
                     ⌊ 
                     
                       
                         170 
                         + 
                         34 
                       
                       2 
                     
                     ⌋ 
                   
                 
               
               } 
             
             . 
           
         
       
     
     S 302 . A three-dimensional trajectory solution set of the target orbit is generated according to the state vector combination corresponding to various time points. 
     S 303 . The 3D trajectory solution set of the target orbit is projected to the instantaneous observation image plane according to the measurement status of the observation platform to obtain a 2D trajectory solution set. 
     Find the corresponding target fitting positions at the same time points and combine them to obtain 2×2×34 pairs of positions and velocity and time combination ({right arrow over (r)} t     s   ′, {right arrow over (v)} t     s   ′, t s ). Calculate the corresponding number of root of orbits {Ô m =(a m , e m , incl m , Ω m , ω m , M m , Ω m ); m=1, 2, 3, . . . , N sv }, the corresponding target orbits {{circumflex over (l)} m =({circumflex over (x)} t     n     m ,ŷ t     n     m ,{circumflex over (z)} t     n     m , t n , n=0, 1, 2, 3 . . . k); m=1, 2, 3, . . . , N sv }, and the trajectory of the target orbit projected onto the image plane: {l*=(U t     n   , V t     n   *, t n , n=0, 1, 2, 3 . . . k); m=1, 2, 3, . . . , N sv }. The obtained target predicted orbit is as shown in  FIG.  7   , wherein the mark indicates an estimated orbit, and mark o indicates a target simulated orbit. 
     S 401 . The derivative error and position error for each two-dimensional trajectory are calculated. 
     The formula for calculating the derivative error of the m-th two-dimensional trajectory is as follows: 
     
       
         
           
             
               Δ 
               ⁢ 
               
                 uv 
                 
                   ′ 
                   ⁢ 
                   m 
                 
               
             
             = 
             
               
                 ∑ 
                 
                   n 
                   = 
                   0 
                 
                 k 
               
                 
               
                 [ 
                 
                   
                     
                       ❘ 
                       &#34;\[LeftBracketingBar]&#34; 
                     
                     
                       
                         ( 
                         
                           
                             2 
                             ⁢ 
                             
                               t 
                               n 
                             
                             × 
                             
                               
                                 U 
                                 ^ 
                               
                               
                                 ″ 
                                 ⁢ 
                                 m 
                               
                             
                           
                           + 
                           
                             
                               U 
                               ^ 
                             
                             
                               ′ 
                               ⁢ 
                               m 
                             
                           
                         
                         ) 
                       
                       - 
                       
                         ( 
                         
                           
                             2 
                             ⁢ 
                             
                               t 
                               n 
                             
                             × 
                             
                               U 
                               
                                 ″ 
                                 * 
                               
                             
                           
                           + 
                           
                             U 
                             
                               ′ 
                               * 
                             
                           
                         
                         ) 
                       
                     
                     
                       ❘ 
                       &#34;\[RightBracketingBar]&#34; 
                     
                   
                   + 
                   
                     
                       ❘ 
                       &#34;\[LeftBracketingBar]&#34; 
                     
                     
                       
                         ( 
                         
                           
                             2 
                             ⁢ 
                             
                               t 
                               n 
                             
                             × 
                             
                               
                                 V 
                                 ^ 
                               
                               
                                 ″ 
                                 ⁢ 
                                 m 
                               
                             
                           
                           + 
                           
                             
                               V 
                               ^ 
                             
                             
                               ′ 
                               ⁢ 
                               m 
                             
                           
                         
                         ) 
                       
                       - 
                       
                         ( 
                         
                           
                             2 
                             ⁢ 
                             
                               t 
                               n 
                             
                             × 
                             
                               V 
                               
                                 ″ 
                                 * 
                               
                             
                           
                           + 
                           
                             V 
                             
                               ′ 
                               * 
                             
                           
                         
                         ) 
                       
                     
                     
                       ❘ 
                       &#34;\[RightBracketingBar]&#34; 
                     
                   
                 
                 ] 
               
             
           
         
       
       
         
           
                                
             
               
                 
                   
                     
                         
                       U 
                     
                     ^ 
                   
                   
                     t 
                     n 
                   
                   m 
                 
                 = 
                 
                   
                     
                       
                         U 
                         ^ 
                       
                       
                         ″ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       
                         U 
                         ^ 
                       
                       
                         ′ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     
                       U 
                       ^ 
                     
                     0 
                     m 
                   
                 
               
               ⁢ 
               
 
                                  
               
                 
                   
                     V 
                     ^ 
                   
                   
                     t 
                     n 
                   
                   m 
                 
                 = 
                 
                   
                     
                       
                         V 
                         ^ 
                       
                       
                         ″ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       
                         V 
                         ^ 
                       
                       
                         ′ 
                         ⁢ 
                         m 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     
                       V 
                       ^ 
                     
                     0 
                     m 
                   
                 
               
               ⁢ 
               
 
                                  
               
                 
                   U 
                   
                     t 
                     n 
                   
                   * 
                 
                 = 
                 
                   
                     
                       U 
                       
                         ″ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       U 
                       
                         ′ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     U 
                     0 
                     * 
                   
                 
               
               ⁢ 
               
 
                                  
               
                 
                   V 
                   
                     t 
                     n 
                   
                   * 
                 
                 = 
                 
                   
                     
                       V 
                       
                         ″ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                       2 
                     
                   
                   + 
                   
                     
                       V 
                       
                         ′ 
                         * 
                       
                     
                     ⁢ 
                     
                       t 
                       n 
                     
                   
                   + 
                   
                     V 
                     0 
                     * 
                   
                 
               
             
           
         
       
     
     wherein U″*, V″*, U′*, V′* are the quadratic polynomial fitting parameters of the two-dimensional trajectory {l*=(U t     n   , V t     n   *, t n , n=0, 1, 2, 3 . . . k)}, Û t     n     m ,{circumflex over (V)} t     n     m  is the corresponding image plane coordinates of the estimated trajectory l m ′ at the time point t n , m represents the m-th two-dimensional trajectory, and k is the total number of observation time points. 
     The camera parameters (orientation) remain the same, and the phase plane coordinates (U t     n   *, V t     n   *) of various time points with respect to the observation of target can be obtained. Quadratic polynomial fitting is performed on the observation data for two times in U and V directions respectively, thereby obtaining the second term and the fourth term. 
     The formula for calculating the position error of the m-th two-dimensional trajectory is as follows: 
     
       
         
           
             
               
                 Rt 
                 m 
               
               = 
               
                 
                   ∑ 
                   
                     n 
                     = 
                     0 
                   
                   k 
                 
                   
                 
                   
                     D 
                     
                       t 
                       n 
                     
                     m 
                   
                   
                     Dc 
                     
                       t 
                       n 
                     
                     m 
                   
                 
               
             
             ⁢ 
             
 
             
               
                 D 
                 
                   t 
                   n 
                 
                 m 
               
               = 
               
                 
                   
                     
                       ( 
                       
                         
                           
                             U 
                             ^ 
                           
                           
                             t 
                             n 
                           
                           m 
                         
                         - 
                         
                           U 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                   + 
                   
                     
                       ( 
                       
                         
                           
                             V 
                             ^ 
                           
                           
                             t 
                             n 
                           
                           m 
                         
                         - 
                         
                           V 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                 
               
             
             ⁢ 
             
 
             
               
                 Dc 
                 
                   t 
                   n 
                 
                 m 
               
               = 
               
                 
                   
                     
                       ( 
                       
                         
                           U 
                           c 
                           m 
                         
                         - 
                         
                           U 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                   + 
                   
                     
                       ( 
                       
                         
                           V 
                           c 
                           m 
                         
                         - 
                         
                           V 
                           
                             t 
                             n 
                           
                           * 
                         
                       
                       ) 
                     
                     2 
                   
                 
               
             
           
         
       
     
     wherein, D t     n     m  is the image plane distance between the predicted value and the observed value, Dc t     n     m  is the image plane distance between the intersection point of the estimated trajectory and the observed trajectory and the observed point, t n  is the observation time, Û t     n     m ,{circumflex over (V)} t     n     m  is the corresponding image plane coordinate of the estimated trajectory l m ′ at the time point t n , U t     n   *V t     n   * is the corresponding image plane coordinate of the observed trajectory l* at the time point t n , and U c   m V c   m  is the image plane coordinate of the intersection point of the estimated trajectory l m ′ and the observed trajectory. 
     S 402 . The derivative error and the position error are taken into comprehensive consideration to select the optimal two-dimensional trajectory, and the number of root of orbit corresponding thereto is used as the optimal estimated value. 
     For the m-th orbit, the derivative error and position error are sorted from small to large, respectively, and the ranking sequence numbers Rank vel,m  and Rank dis,m  are obtained, the two sequence numbers are between 1˜N sv , and the sum of the two is calculated as the final error evaluation.
 
Rank m =Rank vel,m +Rank dis,m  
 
     For the final error evaluation, take the smallest value m opt  as the optimal two-dimensional trajectory: 
     
       
         
           
             
               m 
               opt 
             
             = 
             
               
                 
                   arg 
                   ⁢ 
                      
                   min 
                 
                 m 
               
               ⁢ 
                  
               
                 { 
                 
                   
                     
                       Rank 
                       m 
                     
                     ; 
                     
                       m 
                       = 
                       1 
                     
                   
                   , 
                   2 
                   , 
                   
                     3 
                     ⁢ 
                         
                     … 
                     ⁢ 
                         
                     
                       N 
                       sv 
                     
                   
                 
                 } 
               
             
           
         
       
     
     wherein, N sv  is the number of root of orbits obtained based on the combination of the target position and the velocity data corresponding to various time points. 
     S 403 . The number of root of orbit corresponding to the optimal two-dimensional trajectory is output, and the determination of the initial orbit is completed. 
     As shown in  FIG.  8   , the mark represents all estimated orbits, ▴ represents a preferred estimated orbit, and mark o is a target simulated orbit. 
     The above are only the preferred embodiments of this application, but the scope of this application is not limited thereto. Any changes or replacements that can be easily derived by person skilled in the art based on the technical scope disclosed in the present disclosure should fall within the technical scope disclosed in this application. Therefore, the scope sought to be protected by this application shall be subject to the scope of the claims.